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PREFACE 


This  publication  is  a  revised  and  expanded  edition  of  a  technical  report  written  during  the 
Fall  of  2008.  That  technical  report  was  produced  during  a  flurry  of  activity  hurriedly  conducted  to 
evaluate  a  body  of  technical  work  that  was  relatively  new  to  this  organization.  Out  of  practical 
necessity,  the  research  and  code  development  performed  for  this  project  had  to  be  accomplished 
with  haste  to  meet  the  constraints  of  schedule.  Years  later,  it  seems  that  the  old  report  is  in  need  of 
much  improvement;  also,  a  few  additional  concepts  need  to  be  addressed.  Specifically,  the 
discussion  of  confidence  regions  has  been  clarified  and  made  more  complete.  The  exposition  on 
confidence  ellipses  is  quite  new,  and  in  my  view,  stands  as  an  important  addition  to  the  theory 
behind  sensitivity  testing.  Also,  the  Logit  method  based  upon  the  logistic  distribution  has  been 
derived  in  some  detail.  Additional  examples  have  been  added  to  the  exposition,  and  some  effort 
has  been  invested  in  placing  a  more  substantial  level  of  detail  within  these  exercises.  Hopefully, 
this  information  will  be  of  help  to  the  reader.  Naturally,  as  with  any  first  version  manuscript 
containing  lots  of  equations,  the  errata  list  was  quite  long  enumerating  typographical  mistakes, 
improperly  transcribed  equations,  etc.  A  fair  amount  of  work  has  been  performed  to  correct  these 
mistakes  and  improve  the  quality  of  the  derivations.  I  hope  that  I  have  had  some  success  in  cleaning 
up  the  manuscript,  at  least  in  as  much  as  is  possible  for  someone  who  passed  typing  class  with  a 
score  of  thirty  words  per  minute  with  five  mistakes. 

Statistics  is  a  tough  topic,  so  why  expend  so  much  effort  on  this  report?  (It  is  a  valid  point; 
this  discipline  is  a  tough  study,  and  I  have  spent  a  lot  of  time  on  it).  Well,  it  is  important  as  is 
evidenced  by  its  application  in  the  pharmaceutical  industry,  specifically  in  qualifying  new 
medications.  Of  course,  there  are  many  other  applications  for  these  ideas,  particularly  for  “one- 
shot”  items  that  absolutely  have  to  work  right  when  they  are  used.  Secondly,  this  material  is  very 
interesting  because  it  integrates  different  concepts  in  probability.  This  integral  view  of  probability 
forces  its  practitioner  to  understand  more  of  the  science  behind  the  equations.  In  many  cases, 
people  gravitate  toward  a  deterministic  view  of  the  world  where  the  most  of  one’s  day  is  thought 
of  in  terms  of  absolutes.  That  is  to  say,  one  has  for  instance  a  fixed  amount  of  pocket  change,  37 
minutes  left  to  the  end  of  the  workday,  or  maybe  93  days  left  until  Christmas.  When  you  begin  to 
think  in  terms  of  probability,  a  new  paradigm  must  be  adopted.  Instead  of  thinking  of  these  fixed 
quantities,  we  must  think  instead  of  how  these  quantities  are  distributed  with  respect  to  one  or 
more  variables.  Quantities  are  addressed  in  terms  of  average  values  and  probabilities  instead  of 
exact  values.  This  idea  is  important  when  predicting  the  behavior  of  a  system  with  some  level  of 
certainty.  All  measurements  are  subject  to  uncertainty,  and  these  uncertainties  are  propagated 
through  to  a  system’s  output  with  a  modified  level  of  uncertainty.  This  idea  is  nicely  illustrated  by 
sensitivity  analysis  where  an  uncertainty  governed  by  a  continuous  probability  distribution  is 
integrated  with  the  discrete  binomial  distribution.  The  statistical  properties  of  the  combined 
distribution  are  determined  by  the  algorithms  described  in  this  paper.  The  examples  are  intended 
to  help  the  reader  understand  the  main  concepts. 
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SUMMARY 

This  report  is  expository  in  nature  and  is  intended  to  explain  the  basic  concept  of  a 
sensitivity  experiment  as  well  as  the  analysis  of  its  data.  These  experiments  are  widely  applicable 
to  engineering  problems  that  involve  binary  (pass/fail  of  go/no-go)  outcomes.  We  begin  by 
introducing  the  simple  idea  behind  a  sensitivity  experiment;  then  we  describe  the  basic  “Up  and 
Down”  testing  method.  Particular  emphasis  is  placed  upon  the  parameter  this  procedure  is  intended 
to  identify.  Next,  we  describe  specific  Probit  analyses,  for  analyzing  sensitivity  test  data.  A 
specialized  version  of  this  scheme  is  derived  for  stable  digital  computation.  Confidence  region 
estimation  is  discussed  along  with  an  analysis  of  variance.  The  Logit  method  is  also  presented  as 
an  alternative  to  the  Probit  method.  A  set  of  example  problems  are  solved;  our  results  are  compared 
with  archival  solutions. 
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1  INTRODUCTION 


Statistics  is  a  complicated  science.  Having  terrorized  dormitories  and  classrooms  across 
the  country,  its  late  night,  torturous  study  sessions  have  earned  it  a  dubious  “GPA-busting” 
reputation  along  with  the  forlorn  appellation  Sadistic s  among  captive  students.  Although  it  is  a 
division  of  mathematics,  it  stands  aside  from  the  mainstream  of  this  discipline.  In  fact,  its  more 
skillful  practitioners  refer  to  themselves  as  probabilists  instead  of  mathematicians.  Statistics  can 
also  incite  chagrin  among  practicing  engineers  and  scientists.  In  a  recent  meeting  with  the  author, 
this  assertion  was  reaffirmed  as  one  professional  uttered  the  phrase,  “I  hate  statistics /” 

Why  does  statistics  evoke  such  negative  emotion?  There  is  a  reasonable  answer  for  this 
question.  Although  many  professionals  are  required  to  use  it,  they  lack  the  theoretical  background 
and  experience  required  in  order  to  understand  how  statistics  works.  It  is  a  branch  of  mathematics 
that  differs  from  subjects  like  calculus  or  linear  algebra.  In  these  subjects,  an  individual  can  leam 
a  few  basic  ideas  and  apply  them  in  a  nearly  “cookbook”  way  to  solve  problems.  Statistics  is  not 
so  cooperative.  Statistical  or  stochastic  theory  is  deeply  buried  in  advanced  mathematical  theory 
involving  topics  such  as  Lesbegue  Integration,  Measure  Theory  and  Functional  Analysis. [1]  If 
these  topics  are  well  understood,  then  probabilistic  theory  is  more  accessible,  and  with  practice 
statistics  becomes  understandable.  Unfortunately,  few  people  have  the  time  and  patience  required 
to  study  and  master  these  disciplines.  Hence,  statistical  theory  remains  arcane,  and  its  mastery 
eludes  all  but  its  most  diehard  practitioners. 

The  difficulty  associated  with  mastering  statistical  inference  presents  a  true  dilemma. 
Statistics  is  an  extremely  applied  science.  It  has  an  extensive  number  of  applications  in  practically 
every  field  of  engineering  or  scientific  endeavor.  But  in  order  to  obtain  good  computational  results, 
great  care  is  required  while  converting  abstract  statistical  theories  into  practice.  After  having 
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completed  graduate  study  in  stochastic  processes,  the  author  now  agrees,  at  least  in  part,  with  his 
venerable  instructor  and  thesis  supervisor.  Statistics  is  best  taught  and  understood  through  the  use 
of  clearly  explained  examples  along  with  carefully  planned  excursions  into  probability  theory.  We 
have  attempted  to  take  this  approach  in  the  discussions  that  follow.  He  who  dares  to  venture 
directly  into  the  world  of  stochastic  theory  is  doomed  to  become  lost,  perhaps  forever.  In  other 
words,  when  on  this  safari,  you  need  a  guide. 

This  report  is  designed  to  be  a  guide,  of  sorts.  It  focuses  on  analytical  methods  used  to 
process  binary  (Bernoulli)  and  binomial  statistical  trials.  A  single  binary  trial  or  experiment  has 
only  two  possible  outcomes,  either  a  success  or  a failure.  A  binomial  trial  may  be  thought  of  as  a 
series  of  binary  trials  taken  at  the  same  “level”.  The  result  of  a  binomial  trial  consisting  of  n  binary 
experiments  is  say,  p  successes  and  ( n  -  p)  failures  with  a  crude  probability  of  success  of  p/n  for  a 
particular  binomial  trial. [2]  A  binary  trial  is  a  special  case;  it  is  a  binomial  trial  with  n  =  1  and  p 
can  be  either  0  or  1.  As  it  happens,  Sensitivity  Tests  are  characterized  by  a  mixture  or  series  of 
binary  and  binomial  trials.  For  military  applications,  these  tests  have  a  great  deal  of  utility  since 
many  items  of  military  hardware  may  be  used  only  once  and  may  only  be  evaluated  by  pass/fail 
criteria.  Many  munitions  can  be  thought  of  in  this  way.  For  this  reason,  we  describe  the  basic 
features  of  a  sensitivity  test  in  the  next  section  of  this  report. 

1.1  Sensitivity  Tests 

A  general  sensitivity  test  consists  of  a  finite  series  of  either  binary  or  binomial  experiments 
(or  a  mixture  of  the  two).  Each  experiment  is  performed  at  a  level  defined  by  distinct  value  of  an 
explanatory  variable. [3]  The  value  of  the  explanatory  variable  directly  corresponds  to  the 
measurable  dosage  for  this  level.  The  dosage  is  related  to  the  stimulus,  an  unmeasurable  quantity 
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that  drives  the  outcome  of  the  experiment  (or  trial)  at  the  chosen  level.  [4]  The  stimulus  may  be 
thought  of  as  a  mechanistic  (physical,  biological,  etc.)  process  that  results  from  the  dosage.  Each 
binary  trial  has  only  two  possible  outcomes,  a  success  or  a  failure.  If  a  binomial  trial  (consisting 
of  n  binary  trials)  is  conducted  at  a  given  level,  then  without  any  loss  of  generality,  we  can  say  that 
p  successes  and  (n  -  p)  failures  result  at  this  level.  A  basic  assumption  behind  the  sensitivity  test 
is  that  as  the  dosage  (stimulus)  increases,  then  the  probability  of  success  also  increases.  That  is  to 
say,  at  “lower”  levels  of  stimulus,  we  expect  more  failed  trials.  As  the  stimulus  increases  at  higher 
levels,  we  expect  more  successful  trials,  up  to  the  point  where  practical  all  individual  binary  trials 
are  successes  (or  for  binomial  trials,  p  equals  n).[ 5]  Unfortunately,  due  to  the  pass/fail  nature  of 
each  trial’s  outcome,  we  cannot  precisely  determine  the  dosage  that  will  result  in  a  success.  Instead, 
all  that  we  can  do  is  select  a  dosage  and  then,  via  test,  determine  whether  the  critical  dosage  is 
higher  or  lower  than  the  test  dosage. [6]  That  is  to  say,  if  the  test  is  a  success,  the  critical  dosage  is 
less  than  or  equal  to  the  test  dosage.  Conversely,  if  the  test  is  a  failure,  then  the  critical  dosage  is 
greater  than  the  test  dosage.  To  promote  greater  understanding,  let  us  consider  a  basic  example  of 
sensitivity  testing. 

Suppose  that  a  pharmaceutical  manufacturer  has  developed  a  new  antibiotic  to  fight  a 
particular  strain  of  bacterial  pneumonia.  As  a  normal  part  of  the  certification  procedure,  an 
effective  dosage  must  be  estimated  for  this  drug.  The  critical  dosage  is  defined  as  the  mass  (say, 
in  milligrams)  of  the  drug  that  must  be  administered  to  eradicate  all  invading  bacteria  with  no 
excess  antibiotic  remaining  in  the  body.  As  you  may  imagine,  this  dosage  is  impossible  to  measure. 
Also,  the  outcome  of  the  test  is  determined  by  the  stimulus.  The  stimulus  for  the  problem  requires 
a  detailed  knowledge  at  the  microscopic  scale  of  how  the  antibiotic  interacts  with  each  individual 
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bacterium.  This  information  is  not  known  and  cannot  be  determined,  so  we  must  obtain  an  estimate 
of  the  effective  dosage  by  sensitivity  testing. 

To  estimate  the  effective  antibiotic  dosage,  we  select  a  set  of  dosage  levels  for  testing,  for 
example  {0,  10,  20,  30,  40,  50}  milligrams.  Naturally,  this  set  defines  six  dosage  levels  for  the 
sensitivity  test.  We  then  select  a  number  of  healthy  test  animals  (say,  Rhesus  monkeys)  for  use 
and  decide  how  many  are  to  be  tested.  In  as  much  as  is  possible,  the  test  animals  (or  test  articles) 
should  be  chosen  so  that  they  have  the  same  anatomical  and  physiological  (or  physical) 
characteristics. [7]  If  so  desired,  we  can  formulate  the  entire  sensitivity  test  series  as  a  sequence  of 
binary  trials,  one  test  animal  per  trial.  To  conduct  a  trial  in  this  sequence,  we  choose  a  level  and 
infect  a  test  animal  with  the  pneumonia  bacteria.  Then  we  wait  a  pre-specified  amount  of  time  and 
administer  the  dosage  for  the  chosen  level  to  the  test  animal.  In  the  period  of  time  after 
administering  the  drug,  we  determine  the  response;  the  test  animal  either  lives  (a  success)  or  dies 
(a  failure).  The  response  data  is  collected,  and  we  then  move  onto  the  next  level  and  continue 
testing.  At  the  conclusion  of  the  entire  sequence,  we  have  determined  a  number  of  successes  and 
failures  at  each  dosage  level. 

This  example  brings  an  important  fact  to  light.  As  we  move  between  trials,  we  cannot  reuse 
test  articles,  (e.g.,  test  animals).  Why?  Obviously,  if  a  test  animal  dies  at  the  preceding  trial,  it 
cannot  be  reinfected  and  tested.  If  the  animal  survives  the  preceding  trial,  its  physiology  has  been 
altered  by  the  presence  of  the  antibiotic  and  by  the  action  of  its  immune  system.  The  animal’s 
bloodstream,  lymphatic  system  and  pleural  tissues  are  inundated  with  antibodies  to  the  pneumonia. 
As  a  result,  it  cannot  be  equivalently  reinfected,  so  the  results  of  a  subsequent  test  would  be  tainted. 
Once  a  test  article  has  been  exposed  to  the  stimulus,  it  should  not  be  tested  again.  The  output  of  a 
complete  sensitivity  test  series  may  be  represented  as  a  plot  of  percentage  success  {phi )  versus  the 
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dosage  (or  explanatory  variable)  level.  We  observe  that  the  change  in  the  success  percentage  is 
smaller  for  low  and  high  dosages  while  the  greatest  change  occurs  near  the  mean  dosage.  Hence, 
we  assume  that  the  dosage  response  curve  follows  the  cumulative  normal  distribution  function.  [7] 
When  a  sufficient  amount  of  data  has  been  collected,  this  assumption  is  testable. [7]  With  an 
appropriate  number  of  tests,  we  can  determine  the  “50%  point”  for  the  distribution,  the  level  of 
dosage  that  causes  a  successful  response  for  50%  of  all  like  test  articles.  This  value  of  dosage  is 
the  mean  for  the  distribution. [2]  We  are  also  interested  in  how  success  percentages  behave  for 
dosages  further  away  from  the  mean.  The  property  is  termed  as  dispersion  and  is  measured  by  the 
standard  deviation.  The  primary  focus  of  post-test  analysis  resides  in  determining  these  properties 
for  the  distribution. 

1.2  The  Bruceton  Testing  Method 

The  mean  p  and  standard  deviation  a  must  be  estimated  accurately  if  our  analysis  is  to 
have  real  meaning.  For  this  reason,  the  sensitivity  test  procedure  is  designed  to  concentrate 
measurements  around  the  mean  in  order  to  refine  the  estimate.  The  Bruceton  or  “Up  and  Down” 
test  is  commonly  used  for  this  purpose. [6]  The  first  step  in  conducting  a  Bruceton  test  is  to  decide 
the  range  of  dosages  (minimum  to  maximum)  as  well  as  the  dosage  levels.  We  must  also  decide 
the  total  number  of  individual  binary  trials  to  be  performed  during  the  test  series.  In  practice,  we 
require  a  minimum  of  20  trials,  but  twice  that  number  is  recommended  since  under  theoretical 
limitations,  the  effective  sample  size  is  usually  half  of  the  actual  sample  size.  [6]  Table  1  illustrates 
a  notional  Bruceton  test  series.  The  results  of  the  initial  trial  are  recorded  in  column  T2  while  the 
response  for  the  last  trial  is  recorded  in  column  T21 .  For  the  initial  trial,  we  choose  a  dosage  (level 
2  in  the  example)  that  is  believed  to  be  closest  to  the  actual  mean.  If  the  outcome  of  this  trial  is  a 
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Table  1.  Illustration  of  responses  for  a  Bruceton  test  series.  Series  consists  of  20  binary  trials.  The  notation  "x" 
denotes  a  success  while  "o"  denotes  a  "failure".  Six  dosage  levels  are  numbered  1  to  6  from  low  to  high. 


TRIAL 

T 

T 

1 

T 

2 

T 

3 

T 

4 

T 

5 

T 

6 

T 

7 

T 

8 

T 

9 

T 

10 

T 

11 

T 

12 

T 

13 

T 

14 

T 

15 

T 

16 

T 

17 

T 

18 

T 

19 

T 

20 

Lev  1 

o 

Lev  2 

o 

X 

o 

o 

Lev  3 

0 

X 

o 

0 

X 

0 

Lev  4 

X 

X 

X 

o 

o 

Lev  5 

X 

o 

o 

Lev  6 

X 

Table  2.  Summary  of  the  Bruceton  test  results  extracted  from  Table  1. 


Dosage  Level 

1 

2 

3 

4 

5 

6 

Successes  (x) 

0 

1 

2 

3 

1 

1 

Failures  (o) 

1 

3 

4 

2 

2 

0 

Total  Number  Trials 

1 

4 

6 

5 

3 

1 

success,  we  conduct  the  next  trial  at  the  dosage  just  below  the  initial  level.  On  the  other  hand,  if 
the  first  trial  is  a  failure,  we  conduct  the  second  test  at  the  dosage  level  just  above  the  initial  level. 
The  dosage  level  for  the  third  trial  is  determined  in  the  same  manner  but  is  based  upon  the  outcome 
of  the  second  trial.  Dosage  levels  for  the  remaining  trials  are  determined  in  the  same  way,  but  in 
most  cases,  we  do  not  use  all  of  the  data  collected.  The  reason  for  excluding  some  of  the  data  is 
based  upon  our  desire  to  resolve  the  normal  mean.  As  we  stated  above,  Bruceton  testing 
concentrates  measurements  in  the  vicinity  of  the  mean  dosage.  Since  the  change  in  dosage  level 
between  trials  “opposes”  the  sense  of  the  preceding  response,  the  up  and  down  nature  of  the  test 
series  attempts  to  center  data  collection  at  the  location  of  the  mean.  Unfortunately,  if  our  initial 
guess  for  the  mean  dosage  level  is  poor,  we  automatically  introduce  extraneous  data  into  the  test 
series.  Suppose  that  we  have  guessed  an  initial  value  for  the  dosage  that  is  too  low.  Then  for  two 
or  more  tests,  we  will  steadily  increase  the  dosage  level  and  obtain  failed  responses.  For  example, 
see  levels  2  and  3  in  Table  1.  Let  us  assume  that  we  obtain  a  success  on  the  third  trial  (shown  in 
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column  4).  This  sudden  change  in  the  response  (indicated  by  the  red  block  in  Table  1)  is  termed 
as  a  reversal.  It  follows  that  we  avoid  introducing  extraneous  “start-up”  data  by  analyzing  only 
information  obtained  beginning  with  the  second  trial.  In  this  way,  data  collection  tends  to  remain 
near  the  mean.  The  binomial  results  of  this  example  test  series  are  shown  versus  dosage  level  in 
Table  2. 

1.3  Difficulties  with  Data  Reduction 

Sensitivity  testing  does  have  its  caveats.  In  the  first  place,  the  test  data  must  contain  a  zone 
of  mixed  results.  [8]  That  is  to  say,  the  highest  level  at  which  a  trial  fails  must  exceed  the  lowest 
level  observed  for  a  successful  trial.  If  this  condition  is  not  satisfied,  the  equations  for  estimating 
a  become  inconsistent,  and  zero  is  the  only  value  that  can  be  obtained  for  a. [9]  In  the  example 
shown  in  Table  1,  “2”  is  the  lowest  level  corresponding  to  a  success  while  “5”  is  the  highest  level 
containing  a  failure.  Hence,  for  this  example,  we  obtain  a  zone  of  mixed  results  between  levels  “2” 
and  “5”.  The  Bruceton  test  procedure  is  usually  successful  in  establishing  the  mixed  zone. [8] 

Sensitivity  testing  also  possesses  certain  intrinsic  weaknesses.  Since  it  concentrates  on 
resolving  the  distribution  mean,  it  loses  accuracy  near  the  “tails”  of  the  distribution,  i.e.,  the  regions 
near  0%  and  100%  response. [6]  For  reliability  problems,  we  are  most  interested  in  the  region 
corresponding  to  response  values  exceeding  99%.  As  a  result,  care  must  be  taken  when  sensitivity 
tests  are  conducted  with  this  purpose  in  mind.  Due  diligence  must  be  paid  to  the  structure  of  the 
dosage  levels  and  to  the  number  of  trials.  The  chosen  data  analysis  methodology  is  equally 
important  since  the  choice  of  distribution  can  have  a  significant  effect  in  the  high  probability 
region.  In  certain  cases,  the  logistic  distribution  is  chosen  over  the  normal  distribution  due  to  its 
more  conservative  nature  in  the  “tails”,  i.e.,  the  tails  are  longer. [10]  This  shortcoming  emphasizes 
the  importance  of  correctly  calculating  confidence  intervals  for  system  reliability  probabilities  in 
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the  tail  regions.  Secondly,  the  maximum  likelihood  estimation  (MLE)  procedure  used  to  analyze 
the  test  data  generally  requires  a  large  data  set.  For  small  data  sets,  MLE  loses  “efficiency”;  the 
estimates  of  p  and  a  become  biased  and  acquire  excessive  variance.  [9]  Moreover,  MLE  may  not 
even  be  able  to  analyze  certain  small  data  sets.  A  classic  example  is  that  of  a  test  series  without  a 
mixed  zone  of  results. [9]  Note  that  Dixon  and  Mood  have  published  the  disclaimer:  “Measures  of 
reliability  may  be  very  misleading  if  the  sample  size  is  less  than  forty  or  fifty. ”[6]  To  cope  with 
these  potential  sources  of  error,  careful  test  planning  must  be  combined  with  sound  numerical 
estimation  techniques. 
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2  TECHNICAL  APPROACH 

For  experiments  designed  with  a  fixed  dosage  increment,  one  may  obtain  the  mean  p  and 
standard  deviation  cr  (or  perhaps  other  linear  fitting  parameters)  by  plotting  the  percentage  of 
successes  at  each  dosage  level  on  probability  paper. [7]  If  this  data  correlates  well  with  the  normal 
distribution,  p  and  a  may  be  extracted  graphically  from  the  plot.  Unfortunately,  for  certain  types 
of  sensitivity  tests,  the  dosage  cannot  be  precisely  controlled,  and  the  data  may  not  exactly  conform 
to  the  normal  distribution.  The  former  difficulty  is  commonly  encountered  for  system  reliability 
test  series.  To  achieve  the  best  estimates  of  the  mean  and  standard  deviation  for  these  tests,  we 
usually  apply  generalized  MLE  procedures. [4] 

2.1  Maximum  Likelihood  Equation 

The  Method  of  Maximum  Likelihood  is  a  powerful  estimation  procedure  that  was  developed 
by  R.A.  Fisher  during  the  first  two  decades  of  the  Twentieth  Century.  To  adapt  this  method  for 
analyzing  our  sensitivity  tests,  we  apply  theoretical  concepts  from  probability.  Individual  binary 
trials  performed  in  the  course  of  the  test  series  are  effectively  independent  from  the  standpoint  of 
probability.  The  outcome  (or  response)  of  one  trial  has  no  effect  on  the  outcome  of  any  other  trial. 
If  we  envision  that  the  outcome  of  each  trial  is  represented  by  its  own  random  variable,  then  in 
terms  of  the  dosage  (an  explanatory  variable),  these  random  variables  are  identically  distributed. 
That  is  to  say,  they  have  the  same  distribution  function  and  parameters  ( p  and  a ).  Since  the  trials 
are  independent,  the  probability  of  events  occurring  in  two  trials  is  given  by  the  product  of  their 
individual  probabilities.  Suppose  that  a  total  of  m  binary  trials  are  conducted  at  the  ith  dosage  level. 
Let  the  success  probability  at  this  level  be  given  by  p,  and  the  failure  probability  by 

4,  =  1  -  A  (2.1.1) 
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With  these  assumptions,  the  entire  test  series  can  be  envisioned  as  a  sequence  of  Bernoulli 


(pass/fail)  trials.  At  this  point  in  the  exposition,  we  invoke  Garwood’s  mathematical  expressions 
for  determining  a  lethal  dosage  for  a  drug.  [11]  In  this  lexicon,  a  success  is  equivalent  to  a  death 
resulting  from  a  dosage  while  a  failure  is  equivalent  to  survivor.  Moreover,  consider  that  there  are 
m  -  Si  successes  (deaths)  at  a  given  dosage  level  and  Si  failures  (survivors).  As  an  informational 
aside,  note  that  the  s,  are  random  variables.  The  associated  likelihood  function  can  be  written  as 


i-n 


N  fn^ 


ft:  —S:  S: 


\SiJ 


P  '  '<? 


(2.1.2) 


where  N  is  the  number  of  dosage  levels;  the  term  in  the  parentheses  is  a  combination  of  m  trials 
taken  si  at  a  time.[  11]  The  combination  is  used  to  represent  the  independence  of  order  for 
Bernoulli  trials  defined  at  the  same  level.  Equation  (2.1.2)  is  the  product  of  binomial  probability 
distributions  formulated  at  the  dosage  levels. [2]  Applying  the  natural  logarithm  of  the  likelihood 
function  lends  a  great  deal  of  convenience,  i.e., 


f 


f=i 


L  =  Yj  ln  '  +(ni-si)]n(p,)+siln(cli)\ 


\siJ 


(2.1.3) 


where  L  =  ln(E  ).  By  maximizing  L.  the  natural  log  of  the  likelihood  function,  we  effectively 
maximize  L  .  Equation  (2.1.3)  is  referred  to  as  the  log- likelihood  function. 


2.2  Dosage  Dependent  Probabilities 

To  estimate  the  “mean  dosage”,  the  dosage  defined  at  the  level  where  the  probability  equals 
0.5,  many  probabilistic  concepts  must  be  brought  together.  As  a  result,  it  is  very  easy  to  lose  sight 
of  important  theoretical  details.  One  may  note  that  we  have  done  nothing  to  connect  the  MLE 
expression  to  the  dosage  (or  explanatory  variable).  Equation  (2.1.2)  is  cast  in  the  form  of  a 
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binomial  probability  distribution,  but  one  may  note  that  /?,  and  cp  are  not  yet  related  to  the  dosage 
variable.  These  probabilities  are  related  to  the  dosage  variable  though  the  use  of  a  link  function. [A] 
Several  link  functions  are  available,  but  we  apply  the  normal  probability  distribution  function 
commonly  used  in  the  “probit”  method.  [12]  The  normal  distribution  is  based  upon  a  continuous 
probability  density  function  that  must  be  integrated  with  respect  to  the  explanatory  variable  ( xt  at 
the  ith  level)  in  order  to  compute  p,  where 


Pi  =pM 


clx 


(2.2.1) 


Note  that  by  this  definition,  pi  is  the  probability  that  XE  (-00,  xi )  .  When  (2.2.1)  is  substituted 

into  (2.1.3),  we  obtain  an  expression  for  the  log- likelihood  function  that  explicitly  depends  on  the 
dosage  level.  To  promote  some  simplicity,  we  define  the  variable 


K*)  = 


x-  jU 


o 


(2.2.2) 


and  we  can  rewrite  p,  as 


th 

r 

1  exp 

— 

J— 00 

1  2) 

dt 


(2.2.3) 


This  equation  is  now  in  the  form  of  the  standard  normal  distribution. [9]  Although  (2.1.3),  (2.2.2) 
and  (2.2.3)  are  correct  from  the  standpoint  of  theory,  (2.2.2)  must  be  placed  in  a  different  form  to 
support  the  estimation  procedure. 


2.3  Garwood’s  Method 

In  the  course  of  this  research  project,  the  author  has  encountered  a  number  of  different 
solution  procedures  for  (2.1.3),  (2.2.2)  and  (2.2.3).  There  is  similarity  between  these  methods;  here 
Garwood’s  method  is  presented. [11]  We  begin  this  discussion  with  an  alternate  form  for  t ,  i.e., 
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(2.3.1) 


x  —  ju  _ 

t  = - =  a  +  /?  x 


This  expression  is  just  a  linear  polynomial  in  x  ,  but  it  offers  an  advantage  from  the  standpoint  of 
differential  calculus.  From  (2.3.1),  it  is  easy  to  show  that 


a  1 

u  = - ;  a  =  — 

P  P 


(2.3.2) 


In  the  light  of  (2.1.3)  and  (2.3.1),  critical  points,  local  maxima  and  minima  of  L  ,  may  be  cast  in 
ordered  pairs  (  a,  /? )  ;  they  may  be  determined  by  solving  the  system  of  equations: 


—  =  0;  —  =  0 

da  dj3 


(2.3.3) 


By  differentiating  (2.1.3),  we  may  show  that 


dL  _■ y.  n,  s,  dp, 

dO  i=i  p,q,  \l  n,  J  d6 


(2.3.4) 


where  e  e  (a,j3) .  In  this  case,  we  have  noted  that 


dq,  _  dp, 
80  80 


(2.3.5) 


According  to  (2.2.3),  p,  is  a  function  of  ti ,  so  with  the  use  of  the  chain  rule,  we  may  show  that 


dpL=dtLdpL=dtL  ,  , 
56  86  dt  86PiKi ’ 


(2.3.6) 


By  using  the  process  of  differentiation  with  respect  to  a  parameter  to  (2.2.3),  we  have  that 


1  (  t  I 

P' ^  expl  -J-\  =  Zi;  q,  (t,  )  =  ~z. 


(2.3.7) 


The  derivatives  of  U  ,  in  the  form  of  dtt  led,  arc  calculated  as 


%L  =  \.  d±L  =  x 

5a  ’  d/3 


(2.3. S) 
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Then  (2.3.5),  (2.3.6),  (2.3.7)  and  (2.3.8)  may  be  applied  to  show  that 


dPi  _  .  dpi  _  _  dqi  _  _  3^.  _ 

3a  Z/  ’  dp  x‘Zi  ’  3a  Z/  ’  dp  XiZi 


By  substituting  (2.3.9)  into  (2.3.4),  we  obtain 

dL  N 


r  \ 


_  _  y 

da  tr  /V/, 


<7, 


v  y 


SL  _  y-  _n , 


dp  m  PAi 


\-'-d 

v  *,•  y 


(2.3.9) 


(2.3.10) 


(2.3.11) 


Following  Garwood  [11],  define  (T,  as 


g  =  g(0=^- 


y  _  a 


/V/ 


Z; 


V  'V 


Hence,  (2.3.10)  and  (2.3.11)  become 


3L 

3a  m  ^ 

3L  A 
Sfi~^g‘x‘ 


(2.3.12) 


(2.3.13) 


(2.3.14) 


It  is  important  to  derive  the  second  partial  derivatives  of  the  log-likelihood  function;  (2.3.13)  and 
(2.3.14)  are  helpful  in  this  process.  With  the  use  of  (2.3.12),  we  observe  that 


d2L  3 


N  N 

-I 


*r  ,Ji  _  Y  y 


3a"  3a  3a  cltj  da 


(2.3.15) 


By  using  (2.3.8),  we  obtain 


dL 

da 2 


=If: 


(2.3.16) 


1=1 


In  a  similar  manner,  the  remaining  second  partial  derivatives  are  derived  as  follows. 
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and 


d-L  _  a 

dad/3  da 


d~L  _  d 
d/3  2  ~  d/3 


f  dl? 
v¥y 


r  dO 


f  N  \ 


da 


dt 


a 


v  <=1  y 


6  N  \ 


c)cc  i=\  dt 


(2.3.17) 


i= 1 


v  ’  I  ^  8tj  dgt  2  » 

s/?  “§A'  a/?  ^ 


(2.3.18) 


1=1 


In  order  to  evaluate  (2.3.16)  through  (2.3.18),  an  expression  for  (T;  ;  begin  with  (2.3.12),  and  apply 
the  quotient  rule. 

PAi  ini  Zi  (q,  -  s,  /  >h )}  -  >h  z,  (q,  -  st  /  nt  f  pq, ) 


Si=- 


(pa,)2 


Gi=- 


pa, 


+  Z: 


V  niJ 


<k~- 

\  niJ 


(pa,)2 


'  s  3 

V  niJ 


C pq'i  +  pa, ) 


The  application  of  (2.1.1)  and  (2.3.7)  yields 


Si=- 


PAi 


zf  +  z' 


7, - 

V  nij 


niZi 


(pa,? 


9i-- 

V  niJ 


(-ZiPi+Z#) 


With  careful  algebraic  simplification,  we  obtain 


(  s  ^ 

Q\  = - —  +  7,  ~  ~ 

PAi  PAi  V  n,  j 


p?q 


S; 


i  1  i  V 


"i  J 


PA t 


V  niJ 


For  the  final  three  terms  in  (2.3.22),  we  can  extract  a  common  factor,  i.e., 


Si  = 


7, 


si 

z\ 

zt 

z= 

+  — 

ni; 

Pi 

7,_ 

(2.3.19) 


(2.3.20) 


(2.3.21) 


(2.3.22) 


(2.3.23) 


PAi  PAi 

We  have  values  or  expressions  for  each  of  the  terms  in  (2.3.23)  with  the  exception  of  z[ .  An 
equation  for  this  term  is  easily  obtained  from  (2.3.7).  Observe  that 
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(2.3.24) 


With  the  substitution  of  (2.3.23)  and  (2.3.24)  into  (2.3.16)  through  (2.3.18),  we  arrive  at 
Garwood’s  formulas  for  the  second  partial  derivatives  of  the  log-likelihood  function. 


2.4  Bernoulli  Trial  with  an  Endpoint  Probability 

Garwood’s  method  provides  a  set  of  exact  formulas  for  the  partial  derivatives  of  the  log- 
likelihood  function  with  respect  to  parameters  a  and  y9.[l  1]  We  have  successfully  employed  these 
formulas  to  solve  test  problems,  particularly  Garwood’s  examples.  Unfortunately,  problems  arise 
when  this  method  is  applied  to  sensitivity  tests  that  are  comprised  of  a  mixture  of  binomial  and 

binary  trials.  For  endpoint  probabilities,  those  near  zero  and  unity,  terms  such  as  (T,  and  C,\ 
become  undefined  due  the  presence  of  the  factor  pap  existing  in  the  denominator.  For  endpoint 
probabilities,  Pflj  =  0 ;  this  situation  is  routinely  encountered  for  levels  in  the  explanatory  variable 
characterized  by  a  single  Bernoulli  trial  (or  its  binomial  equivalent).  After  observing  this  difficulty, 
we  elected  to  revisit  the  derivation  and  search  for  an  alternative  form  of  derivatives  dLl  da  , 
dL/dj3 ,  d2L/da2 ,  d2L/dad/3  and  d2L/d/32 .  An  alternative  formulation  is  easier  to  control 
near  the  endpoints. 

This  derivation  begins  by  recalling  (2.1.3)  for  the  log-likelihood  function. 


*■=1 


(2.4.1) 


As  we  differentiate  with  respect  to  a,  we  do  so  term  by  term. 


8L_ 

da 


da  j~i  U 


5  V  M  (  Yi  i 

z.  +  \ni  -  si  /  In  Pi  +  si  In  qt 


(2.4.2) 
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dL_ 

da 


(2.4.3) 


■=I 


1=1 


(«/-■')) 


d  In/;,.  din  q 

'  +  S; 


da 


da 


dL  _  y  | 

da  dor 


fat-s,) 


1  clp .  ^  1  dql 


Pi  dtj  '  q,  dtj 


(2.4.4) 


This  expression  can  be  simplified  by  employing  (2.3.7)  and  (2.3.8)  to  obtain 


dL^y 

da 


Pi  <h 


(2.4.5) 


By  similar  means,  we  can  show  that 


5L  A 

-  =  /  •*; 

s/?  tr 


Pi  Qi 


(2.4.6) 


The  form  of  equations  (2.4.5)  and  (2.4.6)  is  quite  interesting;  the  summand  contains  two  terms, 
and  each  of  these  terms  contains  a  ratio  of  the  probability  density  function  to  its  cumulative 
distribution  function  both  evaluated  at  U  (or  x ,•  under  inverse  transformation).  The  behavior  of  the 
first  (second)  term  ~  Zi  Ipi  is  questionable  as  pi  — >  0  .  Conversely,  the  behavior  of  the  second  term 
(~  Zi  Iqi )  requires  investigation  as  p,  — >  1  (or  q;  — >  0).  If  we  can  assert  that  z.i  Ip,  tends  to  zero  as 

Pi  ~ >  0 ,  we  can  exclude  this  term  from  (2.4.5)  in  the  limit.  This  issue  is  addressed  by  the  following 
proposition. 

Proposition:  Let  z,  be  the  normal  probability  density  function,  and  let  p  be  its  cumulative 
distribution  function  both  defined  in  the  independent  variable  t.  The  ratio  z/p  is  well  defined  and 
tends  to  zero  as  p  — >  0. 

Proof:  By  examining  equations  (2.3.7)  and  (2.2.3),  z  and  p  jointly  approach  zero  if  and 
only  if  t  — >  -oo.  The  ratio  z/p  may  be  written  as  follows. 
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)(-  r  /  2) 


-W=— expv 

Jexp(-  f2  /2)<7? 


(P-1) 


The  Taylor  series  expansion  for  the  exponential  function, 


00  xl 


exp(x)  =  1  +  x  +  — — l —  =  '^J~ 

2  ,=o  i  ! 


is  a  convergent  power  series  for  x  <  0.  Hence,  it  follows  that 


(P-2) 


exp 


V  2  y 


2  4 

=  1 - —  + 


2-1!  4-2!  8-3! 


+  ...  =  £si+o(t2n) 


(P-3) 


where 


(-1 V 

2'i\ 


(P-4) 


is  also  a  convergent  power  series  for  real  valued  t.  To  evaluate  p,  as  occurs  in  (P-1),  we  compute 
the  anti-derivative  of  (P-3)  term  by  term,  i.e., 


ct 

f  t2) 

1  exp 

J— CO 

t 3  t5 

dt  =t - + 


6-1!  20-2!  56  •  3! 


+  --  =  tfjSi+o(t2n+1)  (P-5) 


where 


S_t#l 

'  2r(2j  +  l)i! 

The  convergence  of  this  series  is  confirmed  by  applying  the  ratio  test.  Observe  that 


(P-6) 


sM 

(-1  r(<T 

2'(2/  +  l)/! 

1 

To' 

+ 

to 

S, 

2'+1  [2(7  +  !)+!](/  +  1)! 

(-Ot) 

2  i  +  5/  +  3 

(P-7) 


From  this  ratio,  we  conclude  that 


r  K 

lim 


=  0 


(P-8) 
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Hence,  the  power  series  (P-5)  converges  for  any  value  of  t.  Since  both  of  the  power  series 


representations  for  the  functions  in  (P-1)  are  convergent,  both  have  finite  limits,  say  f(t)  for  the 
numerator  and  g  (/ )  for  the  denominator,  so  for  any  real  t 

£  -  fill  XM  ,  '  _  fill  (P-9) 

P  tJStg{t)-7giJ)  tg(t ) 


having  invoked  the  properties  of  the  cumulative  distribution  function  as  t  — >  — °o .  Since  /(/  )  and 
g(t)  are  the  finite  limits  of  series  that  converge  with  the  same  order,  we  have  that 

lim-  =  lim/rl  =  0-  ° 

p-> 0  P  p— >0  t  g[t ) 

By  using  a  proof  similar  to  that  shown  in  the  preceding  proposition,  the  ratio  zJq  also 
approaches  zero  as  p  — >  1  (q  — >  0  ).  It  follows  that  (2.4.5)  and  (2.4.6)  can  be  expressed  as 


N 


U= 1 


and 


dL_ 

dp 


N 

-  Z  x<  s.  z>  Pi=0 

i=l 

N 

Yjxi(ni  ~si)zi  cp  = 0 

.  i= 1 


(2.4.7) 


(2.4.8) 


By  using  the  above  equations  for  endpoint  probabilities,  dL/  da  and  8LI  dp  ,  derivatives  of  the 
log-likelihood  function  remain  well  defined  for  all  possible  values  of  pi  and  qu 

In  addition  to  the  Jacobian  matrix,  we  must  also  use  the  matrix  of  second  partial  derivatives 
or  Hessian  matrix.  The  reason  for  its  necessity  will  be  made  clear  in  the  next  section.  We  would 
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like  to  construct  the  second  partial  derivatives  in  terms  of  the  ratios  Z,  /  p,  and  Zt  /  Cj, .  To  derive 
d2L! da2 ,  we  differentiate  (2.4.5),  i.e., 


v  (  )  5  U  8  z< 

L  \ni~sih - - 


da\pi  J  dayqi 


(2.4.9) 


With  some  careful  mathematics  and  the  use  of  (2.3.7)  and  (2.3.24),  we  can  show  that 


x  f  \  j  f  \  \ 

_d_  ^dt^d_  Zi_  =_t  _  Zi_ 

do  \  Pi  j  dadti{pi)  '{pj  [pi 


(2.4.10) 


a  f  \  a*  j  I  \  \ 

_d_  Zi_  =  dt^d_  =_t  + 

do{qi)  SadtXqJ  '{qj  [qi 


(2.4.11) 


After  substituting  (2.4.10)  and  (2.4.11)  into  (2.4.9),  we  obtain 


UA  Pi  UA  CU 


(2.4.12) 


Note  that  (2.4. 12)  is  characterized  by  the  presence  of  the  ratios  Z,  I  p{  and  Z,  /  q,  but  no  other  terms 
containing  Zi  and  pu  The  mixed  partial  derivative  is  developed  as  follows. 


dad/3  tP  PUJ  '  dfiW 


(2.4.13) 


By  carefully  differentiating  and  using  (2.3.7)  and  (2.3.8),  we  can  show  that 


d  f  zt  |  _  dtj  d  f  zt  |  _  [  zt  j  f  Zi 

dP  i  Pi  J  dp  dti  IpJ  *'*'[/>,■  J  [Pi 


(2.4.14) 


d  Zj  _  dti  d  z{  _  t  Zi  _  Zi 

dpyqj  dpdtXqi)  '  \  \q,  J  U 


(2.4.15) 


By  substituting  (2.4.14)  and  (2.4.15)  into  (2.4.13),  the  mixed  partial  derivative  is  revealed  as 
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d2L 

dad/3 


-2>. 


i=l 


Pi) 


‘■+-, 

Pi) 


~  S; 


.  3,) 


(2.4.16) 


Note  that  the  ratios  found  in  (2.4.12)  recur.  If  (2.4.6)  is  differentiated,  the  remaining  second  order 
partial  derivative  may  be  calculated. 

d2L  d 


d/32  d/3 

By  simplifying,  we  have  that 


Pi  3, 


z, 

-S;—} 


1=1 


(2.4.17) 


d2L  A 

=  Lxi 


dj3z 


i=i 


V 


d/3 


PiJ 


d 


S; 


'  d/3 


f  \ 

h. 

Uv 


(2.4.18) 


Substitution  of  (2.4.14)  and  (2.4.15),  produces 


d2L  ^ 

Lxt 


d/3 


1=1 


' 

f 

(  \ 

(  \ 

2 \ 

' 

( 

(  \ 

(  \ 

2 N 

(  \ 

zi 

Zi 

Zi 

Zi 

\n,  -sj 

h 

+ 

>-sr 

~xi 

h 

— 

> 

V 

{ Pi) 

{ Pi) 

J 

V 

yVi  y 

\Vi  j 

J 

(2.4.19) 


d2L 

d/32 


-Iv 


i=i 


t;+  — 


~  S; 


V' 


Vi) 


h-1 

V  Vi  J 


(2.4.20) 


Pi  A  P<  J 

To  recap,  this  section  presents  an  alternative  formulation  of  Garwood’s  method  that  is  useful  in 
evaluating  endpoint  probabilities.  At  these  endpoints,  the  second  order  partial  derivatives  may  be 
computed  as  follows. 


d2L 
da 2 


d2L 

dad/3 


Pi  = 0 

1=1 

N 

~Tj(ni~si)zi(ti+zi)  Vi  =0 
1=1 

N 

Pi=  0 

1=1 

N 

~Tuxi(ni  ~si)zi(fi  +  zi)  Vi  =0 


(2.4.21) 


(2.4.22) 
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Pi  =  0 


(2.4.23) 


Yjxrhziih  ~  zi) 

i= 1 
N 

-Yjx?(ni~si)zi(ti  +  zi)  v, =0 
1=1 

For  the  endpoint  probabilities  ( p  =  l,  q  =  o )  and  ( p  =  o,  q  =  1 ),  equations  (2.4.7),  (2.4.8),  (2.4.21), 
(2.4.22)  and  (2.4.23)  are  instrumental  in  supporting  the  procedure  for  estimating  //  and  a. 

2.5  Newton’s  Solution  Method 

In  Sections  2.3  and  2.4,  it  was  demonstrated  that  we  may  fit  sensitivity  test  data  to  the 
normal  distribution  through  Fisher’s  Method  of  Maximum  Likelihood. [13]  Moreover,  the 
procedure  requires  that  we  estimate  two  parameters,  the  distribution  mean  and  standard  deviation, 
p  and  er,  respectively.  To  do  so,  we  must  determine  a  and  fi  satisfying  (2.3.3).  Estimates  of  p  and 
a  are  obtained  by  solving  equations  (2.3.3)  cast  either  in  Garwood’s  form  (Section  2.3)  or  an 
alternative  form  (Section  2.4).  The  fitting  procedure  takes  the  form  of  an  iterative  scheme;  rewrite 
a  and  /I  as  follows. 

a  =  a0  +  Aa  ;  /?  =  /?„  +  A/?  (2.5.1) 

Substitute  (2.5.1)  into  (2.3.3);  when  a  Taylor  series  is  constructed  in  variables  a  and  /?,  observe 
that 

|^(a,/9)=|^(a„.,»o)+Aa|%(a0,/?(J)+AJsT2-(a0,/?0)+o(Aa!,AaA/?)  (2.5.2) 
da  da  da~  da  dp 

dl  dT  d~T  d~T  i  \ 

(«o.A)+Aa— -(«„A)+4/»- T(a0,/?0)+OAoA/?,A/;!  (2,5.3) 

dp  dp  da  dp  dp 

An  iterative  estimation  scheme  may  be  retaining  only  the  first  order  terms  in  (2.5.2)  and  (2.5.3). 
In  matrix  form,  we  have  that 
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"  d2L(  ^  82L  (  a^\\*a 

aa2(«°.A)  QaQpbo'h) 

/^(«o.A)  0(«o.A)  - 

5a  dp  op  a  B 


5a 

f^o>A>) 


(2.5.4) 


The  2x2  matrix  on  the  left  side  of  (2.5.4)  is  H ,  the  Hessian  matrix,  and  elements  of  the  Jacobian 
matrix  (a  vector  J  in  this  case)  resides  on  the  right  side.  The  pair  {cc[V  fi{])  is  designated  as  a 

starting  guess;  then  the  pair  (Aa,A/3)  is  an  increment  that  when  added  to (<Z0,/?0)  ostensibly 
forms  an  improved  estimate,  a  new  iterative  value  for  (a,j3),  and  by  transformation  (2.3.2),  a  new 
estimate  for  (p,o).  The  system  (2.5.4)  may  be  solved  for  the  increments  by  inverting  H  and 


multiplying  from  the  left.  The  inverse  2x2  matrix  H  1  is  easily  obtained  as  long  as 

|H(a0 ,  /?o )  I  ^  0 


(2.5.5) 


The  updated  increment  (a A/3)  is  calculated  from  the  equation 


[Aa,A/?f  =  H“'(a0,/?0)*  7(a0,/?0) 


(2.5.6) 


It  follows  that 


[ct,p]  =[a0,/?0f  +H“1(a0,^0)*/(a0,/?0) 


(2.5.7) 


where  [a,/?]7  is  a  vector  that  contains  the  refined  parameter  estimates.  As  with  other  Newton 


approximation  methods,  (2.5.7)  is  easily  structured  as  an  iterative  scheme.  After  a  number  of 
iterations,  when  this  scheme  has  converged  to  the  desired  level  of  accuracy,  the  distribution  mean 
and  standard  deviation  can  be  computed  from  (2.3.2). 
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2.6  Minor  Variant  of  the  Golub-Grubbs  Method 

Garwood’s  method  for  probit  analysis  is  highlighted  in  the  preceding  sections.  This  method 
is  based  upon  a  linear  fit  (2.3.1)  of  the  probit.  Recall  that  the  probit  is  the  argument  to  the  normal 
probability  integral  [8] 

(x-fi)  Y  l  (  r) 

rh(^expN  <2'6'1) 

the  probit  being  defined  as 

t  =  (2.6.2) 

cr 

The  maximum  likelihood  method  can  be  applied  with  the  form  (2.6.2).  For  this  variant  of  the 


method,  let  p{  be  the  success  probability  while  St  is  the  number  of  successes  out  of  trials.  In 
this  case,  the  log-likelihood  function  is  written  as 


*=Z 


+  si\n(pi)+(ni-si)\n(qi) 


(2.6.3) 


To  determine  the  value  of  the  explanatory  variable  x  for  which  the  log-likelihood  function  is 
maximized,  it  is  necessary  to  differentiate  (2.6.3)  with  respect  to  q  and  cr,  i.e., 


dh_ 

dju 


z  *  4-  (in  Pt ) + ( ni  -  si  )4~  (in  % ) 

,•=1  dM  dp 


(2.6.4) 


dL_ 

da 


Y"  ^  Pi)  +  ini 

da 


Si)—  (ln^,) 

da 


(2.6.5) 


Applying  the  chain  rule  along  with  the  derivative  of  the  logarithm, 

dL  s‘  8p’  |  (ni~si)dPi 
dp  m  Pi  dp  q,  da 


(2.6.6) 
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dL  -  V  -s',  dP,  |  (ni-s,  )dPi 


(2.6.7) 


da  i=l  pi  da  q{  da 

To  evaluate  (2.6.4)  and  (2.6.5),  the  following  additional  derivatives  are  required.  From  (2.6.2),  we 
have  that 


dt;  d  (  Xj  -  p  ^ 


dju  dp 


V  cr  J 


1 

a 


(2.6.8) 


dt.  i 


da  dp 


V  cr  J 


c2J 


A 

a 


(2.6.9) 


Additional  derivatives  are  shown  in  (2.3.7)  and  (2.3.24);  therefore, 

dpL=dpLdtL  =  _z±  '  dp_  =Zi_ 
dp  dt ,  dp  a  dp  a 


(2.6.10) 


0pL=  dp^dt^  =  _zp_  _  dq^=zp_ 
da  dt ,  5(7  (7  ’  5cr  cr 


(2.6.11) 


Since  qt  -  1  -  pt , 


Al  =  A.  %=Al 

dp  a  '  dp  a 


(2.6.12) 


By  substituting  these  derivatives  into  (2.6.6)  and  (2.6.7), 


SL  _ 

5//  <7  ,=1 

5L  _  1  ^ 
<9cr  a  f 


+(«,• 

-*,)$- 

V  Pi 

dt) 

-  v,-  — +(«, 

-s,)tA 

Pi 

di 

(2.6.13) 


(2.6.14) 


A  maximum  occurs  in  the  log-likelihood  function  for  the  pair  (// ,  cr)  when 


dL  _  5L  _0 
5//  5cr 


(2.6.15) 
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d2L  d2L  d2L 

It  is  also  necessary  to  determine  the  second  partial  derivatives  — j , - and  — — j .  Based  upon 

dp "  dp  da  da' 

the  ratios  in  (2.6.13)  and  (2.6.14),  the  following  ordinary  derivatives  are  needed. 


d_ 

dt- 


/  \ 

f  \ 

r 

( 

f  \ 

r  \ 

2\ 

zt 

= 

zt 

—  — 

6 

Zi 

+ 

Zi 

UJ 

[Pi) 

V 

{Pi) 

) 

(2.6.16) 


r  \ 

(  X 

r 

( 

(  X 

r  \ 

2\ 

h. 

_ 

A 

_  _ 

t 

A 

A, 

a) 

1 1 

V 

Ah, 

Ah) 

y 

(2.6.17) 


_d 

dt, 


i  \ 


Pi  J 


A 

V  Pi  J 


i~i  |  _  Zj  _ j.2 


f  7  \  f  -  3  2 
'  '~t. 


Pi 


Pi J 


Z, 

\PiJ 


(2.6.18) 


A 

dt 


rt:z ' 


i  V 


CU  J 


A 

v  9i  J 


i  |  _  +  _ p 


f  \  f  \2 
+  61  Zi 


z, 

CU) 


ch) 


(2.6.19) 


Ch  v 

Notice  that  the  ratios  Zi/pi  and  z.i  Iqi  repeatedly  occur  in  these  derivatives.  With  all  of  these 
derivatives  in  mind,  the  second  order  partial  derivatives  may  be  derived  as  follows. 


d2L  ]_ A 

c  ;=1 


dp2 


0  1 

{a 

i 

G 

^1 

l  Pu 

+  {ni~si) 


_d_ 

dp 


dp~ 


s2l  _  1  V 

O'  M 


f  X 

t 

f  X 

t 

1  Zi 

(  'J  Zi 

si 

— 

-{p-SiJ  — 

)pi) 

UJ 

f  \ 

A 

A  J 

06 

dp 


By  substituting  (2.6.8),  (2.6.16)  and  (2.6.17), 


d'L 
dp 2 


1  N 


i= 1 


-  S’,- 


v  6,  y 


/  A 
Z; 


\PiJ 


+ 


Zj_ 

\di  J 


r  \ 
h. 

\di  j 


2  6 


(2.6.20) 


(2.6.21) 


(2.6.22) 


Similarly,  the  mixed  partial  derivative  can  be  developed  as  follows.  Differentiate  (2.6.14)  with 
respect  to  p. 
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d2L  1  A 

a  1 

(VT 

d//0cr  cr 

i 

c. 

#1 

l  Ay 

+ 


(ni~si) 


_d_ 

dju 


( ,  A 


v  A  y 


(2.6.23) 


Hence, 


d2L 


1  N 


-  5, 


djuda  \  P;  ) 

and  by  applying  (2.6.8), 

d2L  _  1  ^ 
djuda  o2 

Substituting  (2.6.18)  and  (2.6.19),  we  obtain 


+  (ni~si) 


ft;Z  " 


V  A  J 


dh_ 

d/u 


(2.6.24) 


/  \ 

t 

f  \ 

6-z,. 

1  /  \l  6Z; 

•3 

1 

7s 

1 

_&5 

l  A  J 

1  l  A  J 

(2.6.25) 


8~L 


i  N 

=Ai 


a//  So-  o"  tr 


5,. 


'Z 


vAy 


-t; 


'z' 


v  A  j 


-t, 


2  3 


Aj 


(a  “  A 


v 

Ay 


y  \ 
Z/_ 

vAy 


+  t  ■ 


r  32^ 
z 


v  y 


(2.6.26) 


The  final  partial  derivative  requires  a  more  tedious  mathematical  development.  By  differentiating 
(2.2.7),  we  obtain 
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Application  of  the  product  rule  produces 
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Simplifying  with  the  use  of  (2.6.9),  we  obtain 
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The  components  of  this  system  a  common  factor  of  1/cr.  Since  the  absolute  value  of  a  can  be  quite 
small,  it  is  numerically  advantageous  to  multiply  the  matrix  and  vector  in  (2.6.31)  by  a2.  Doing 
so  has  the  effect  of  increasing  the  magnitude  of  the  matrix  components,  particularly  the  diagonal 
entries.  Accomplishing  this  change  enhances  the  convergence  of  the  iterative  scheme. 
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2.7  Analysis  of  Confidence  in  the  Estimates 

One  tool  for  analyzing  the  level  of  believability  associated  with  the  estimate  of  an  unknown 
parameter  is  the  confidence  region.  [14]  The  confidence  region  represents  one  of  the  more  arcane 
concepts  in  statistics,  so  it  is  best  illustrated  for  a  single  explanatory  variable  say,  the  “ideal” 
dosage  for  an  antibiotic.  When  the  ideal  dosage  is  administered,  it  should  manifest  a  concentration 
in  living  tissue  that  is  sufficient  for  the  eradication  of  the  pathogenic  bacteria.  Given  the  level  of 
uncertainty  in  dosage  measurement,  administration  and  effectiveness,  we  would  like  to  estimate 
an  effective  range  for  dosage  quantity  and  then  to  attach  to  it  a  probabilistic  level  of  assurance. 
The  dosage  is  a  quantity  represented  in  one  dimension,  so  in  terms  of  the  success  probability,  its 
range  is  designated  as  a  confidence  interval,  and  the  mathematical  process  used  to  determine  the 
interval  is  known  as  an  interval  estimate.  The  end  points  of  the  interval  act  as  random  variables 
while  the  statistical  parameter  theoretically  contained  within  the  interval  is  regarded  as  fixed,  but 
unknown.  [13]  The  concept  of  the  confidence  interval  can  be  confusing,  so  it  is  instructive  to 
explore  why  this  concept  can  present  difficulties. 

Consider  a  physical  or  biological  process  that  exhibits  some  randomness  or  “noise”. 
Suppose  that  the  behavior  of  this  process  can  be  characterized,  in  part,  by  an  unknown,  yet  fixed 
non-ranclom  parameter  9.  A  good  example  of  this  parameter  is  the  mean  of  a  statistical  distribution. 
The  parameter  exists,  but  its  exact  value  cannot  be  known.  Instead,  9  must  be  estimated  via 
experimentation,  perhaps  by  using  a  series  of  independent  experimental  trials.  To  differentiate  it 
from  the  unknown  exact  value,  the  estimated  value  of  9  is  denoted  as  9 .  Based  also  on  information 
provided  by  the  experiments,  an  interval  estimate  \o,  ,  0H  ]  associated  with  0  may  be  calculated. 
For  this  interval,  the  investigator  would  compute  a  probability  p  for  the  interval  and  then  make  the 
following  assertion. 
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CLAIM:  With  lOOp  %  certainty,  the  exact  value  of  the  parameter  9  is  contained  within  \oL  ,0H\. 

It  is  in  consideration  of  this  assertion  that  the  difficulty  arises.  Although  it  sounds 
reasonable,  this  claim  is  false.  Why?  Procedurally,  the  interval  endpoints  (values  of  random 
variables)  are  calculated  based  upon  the  value  of  0.  Thus,  at  the  conclusion  of  the  experiments  (a 
series  of  Bernoulli  of  binomial  trials),  0  is  computed;  then  ()L  and  0H  are  calculated  directly.  As 
they  pertain  to  a  single  set  of  experiments  9 ,  0L  and  (JH  are  fixed.  Since  they  have  been  calculated, 
these  numbers  have  no  random  behavior.  Therefore,  we  cannot  assign  a  probability  to  any  of  these 
values,  so  the  claim  made  above  cannot  be  true.  For  that  reason,  it  is  important  to  understand  the 
true  probabilistic  interpretation  for  this  situation. 

Recall  that  in  the  mathematics  of  probability,  the  interval  endpoints  are  represented  by 
random  variables.  Denote  these  random  variables  as  © L  and  ©  R  .  An  initial  test  provides  a  first 
estimate  for  these  values  (and  for  §).  If  we  repeat  the  same  experiment,  that  is,  the  series  of  trials, 
a  number  of  times,  we  obtain  a  sequence  of  values  for  0,  and  for  &R  .  For  any  one  of  these 
tests,  there  is  nothing  random  left  to  constitute  a  random  interval  endpoint.  The  probabilistic 
aspects  of  the  interval  emerge  only  when  we  view  the  series  of  repeated  tests  in  a  larger  sense. 
That  is,  each  test  becomes  a  separate  realization.  In  this  sense,  the  interval  endpoints  fluctuate 
randomly.  If  100  tests  are  performed,  then  we  obtain  100  random  intervals  [©L  ,  (-) R  ] .  Therefore, 
randomness  is  effectively  restored  to  the  notion  of  the  interval.  With  this  interpretation,  one  may 
wonder  as  to  how  the  faulty  claim  made  above  might  be  replaced  by  a  statistically  sound  assertion. 
Suppose  that  we  are  interested  in  the  parameter  9,  and  consider  a  95%  confidence  interval 
associated  with  9.  The  statistically  sound  assertion  is  made  as  follows. 
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PROPER  CLAIM:  Suppose  that  an  unknown  statistical  parameter  6  is  estimated  in  100  separate, 
repeatable  experiments.  If  95%  confidence  interval  endpoints  are  calculated  for  each  experiment, 
then  the  true  value  of  the  parameter  6  will  be  contained  in  95  of  the  100  intervals.  The  probability 
that  the  true  value  of  0  lies  within  the  interval  computed  for  an  individual  test  is  either  one  or  zero. 

In  the  paragraphs  that  follow,  the  mathematical  development  of  a  confidence  interval  is 
presented.  This  exposition  is  not  rigorous,  but  the  practical  mathematics  behind  interval 
construction  is  emphasized.  Elementary  estimation  theory  is  drawn  upon  as  a  core  resource. 

2.7.1  Success  Ratio  Confidence  Interval 

For  the  success  ratio  (or  success  probability)  confidence  interval,  the  success  probability 
p  is  estimated.  The  associated  interval  [  Pl  ,  Pr\  is  determined  with  a  confidence  level  of  C  where 

C  =  100(l-a)%,  0  <  a  <  1 

If  100  separate,  but  identically  constructed  tests  are  performed,  then  100  estimates  of  p  are 
computed  along  with  100  intervals  [Pl  ,  Pk  \ .  Hence,  by  invoking  the  proper  claim  made  in  the 
previous  section,  100(1-  a)  of  these  intervals  will  contain  the  true  success  probability.  The  interval 
endpoints  may  be  derived  by  using  theory  associated  with  the  binomial  distribution. [15]  The 
starting  point  for  this  derivation  is  the  DeMoivre-Laplace  Theorem  where  it  is  assumed  that  n,  the 
number  of  binomial  trials  conducted  for  the  experiment,  is  large.  Suppose  that  Y  is  the  number  of 
successful  trials  and  that  p  is  the  unknown  success  probability;  then  the  DeMoivre-Laplace 
theorem  states  that  the  ratio  z  defined  by 
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Figure  1  Probability  Distribution  Function  for  the  Standard  Normal  Distribution 
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tends  toward  the  standard  normal  distribution  /V(0,1),  i.e., 
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In  terms  of  probability,  this  theorem  implies  that 


Zal 2  < 


IPd-p) 


<  z 


all 


=  1  -a 


(2.7. 1.1) 


(2.7. 1.2) 


(2.7. 1.3) 


The  similarity  of  (2. 7. 1.1)  to  (2.7. 1.2)  can  be  illustrated  by  realizing  that  for  the  binomial 
distribution,  //  =  np  and  =  n p{\~  p) \  15]  With  these  substitutions,  (2.7. 1.1)  can  be  rewritten 


as 
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Y-p 


(2.7. 1.4) 
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This  expression  is  identical  to  (2.2.2),  the  argument  of  the  standard  normal  distribution.  Denote 
p  =  Y/n  as  the  estimate  of  the  success  probability;  then  (2.7. 1.3)  becomes 


MO,  1 )  is  a  symmetric  function  about  its  zero  mean,  so  a  typical  confidence  interval 
(-  Za/2 , Za/2)  is  also  symmetric  about  the  mean.  The  probability  associated  with  this  interval  is 
obtained  by  integrating  the  density  function  (2. 6. 1.2)  over  (-  Za/2  >  Za/2) .  As  is  illustrated  in  Figure 

1,  the  regions  00  Zan)  and  {zai2 , 00 )  form  the  “tails”  of  the  distribution.  The  tails  are 
important  when  determining  the  probability  or  confidence  coefficient  associated  with  the  interval. 
As  is  shown  in  Figure  1,  the  probability  (area  under  the  curve)  for  the  tail  regions  is  a.  Then  the 
confidence  coefficient  for  the  interval  is  defined  as  100(l-a)%.  For  example,  consider  a  confidence 
coefficient  of  95%;  the  attendant  value  of  a  is  0.05.  It  follows  that  the  probability  in  each  tail  region 
is  all  =  0.025.  The  left  interval  endpoint  on  the  explanatory  z  axis  for  M0,1)  can  now  be 

determined  by  solving  for  -  Za/2  in  the  cumulative  distribution  function  equation 

-v1—  Z\2  expf-  —1  dz  =  —  (2.7.1.6) 

v/2 n  U  y  2  )  2 

Taking  (2.6. 1.5)  into  careful  consideration,/)  is  identified  as  the  unknown  probability  while  p  is 
its  estimate.  By  algebraically  solving  for  p  in  (2.7. 1.3),  an  interval  estimate  for  p  can  be  obtained 
consistent  with  the  interval  probability  l  -  a.  Collett  [4]  derives  two  expressions  for  the  interval; 
the  first  is  derived  as  follows.  Begin  by  rewriting  (2.6. 1.5). 
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(2.7. 1.7) 


The  bounds  of  the  interval  are  identified  by  selecting  the  equality  in  (2.7. 1.7);  by  rearranging  terms, 
we  obtain 
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Additional  algebraic  manipulation  produces  the  quadratic  equation 

+  ;32  =  0  (2.7. 1.9) 

From  an  analysis  of  its  form,  (2. 7. 1.9)  defines  a  parabola. [15]  The  roots  of  this  equation  delimit 
the  endpoints  of  the  confidence  interval.  Formulas  for  the  endpoints  are  calculated  by  using  the 
quadratic  formula.  After  some  simplification,  we  obtain 
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(2.7.1.10) 
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This  formula  contains  the  so-called  standard  error,  se(  p  )  for  the  binomial  distribution  [4],  i.e., 

p( 1  -  P) 


se 


(p)  = 


(2.7.1.11) 


With  the  use  of  (2.7.1.10)  and  (2.7.1.11),  the  interval  endpoints  pu>w  and pmgh  are  written  as 
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(2.7.1.13) 
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These  formulas  deserve  added  commentary.  Recall  that  the  confidence  coefficient  is  based  upon 
M0,1)  and  the  symmetric  interval  (~  Za/2  »^a/2)-  The  interval  delineated  by  (2.7.1.12)  and 
(2.7.1.13)  is  symmetric  about  the  quantity 
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(2.7.1.14) 


where  p  can  be  thought  of  as  an  observed  and  rather  crude  success  ratio  probability.  It  is 
interesting  to  note  that  p  does  not  lie  at  the  center  of  the  interval.  Not  only  has  this  probability 
been  translated  off  of  the  estimate,  it  has  also  been  scaled.  Recall  that  the  DeMoivre-Laplace 
theorem  states  that  the  binomial  distribution  tends  toward  the  normal  distribution  when  n  becomes 
large.  Equation  (2.7.1.14)  is  consistent  with  the  conditions  of  this  theorem.  Note  that  p 
approaches  p  as  n  becomes  very  large.  Equations  (2.7.1.12)  through  (2.7.1.14)  delineate  a  first 
model  (denoted  as  Model  1)  for  the  endpoints  of  the  confidence  interval  for  the  success  probability. 

The  second  model  for  the  success  probability  is  derived  by  recalling  (2.7. 1.7).  In  this  case, 
the  right  hand  side  of  the  equation  is  evaluated  by  substituting  p  ,  the  estimate  for  the  success 
probability.  By  rewriting  this  expression  as  an  equality  at  the  endpoint,  we  obtain 
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(2.7.1.15) 


Hence,  with  the  use  of  (2.7.1.11),  we  obtain  the  second  model  (denoted  as  Model  2),  i.e., 


Plow  ~  P~  Zal2  Se  (P) 


(2.7.1.16) 
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Pmgh  =P-  zal2  se  (p) 


(2.7.1.17) 


Model  2  is  more  commonly  used  for  delineating  the  endpoints  of  this  confidence  interval. 

2.7.2  Success  Probability  Confidence  Interval  in  the  Explanatory  Variable 

For  sensitivity  testing,  the  binomial  distribution  remains  a  critical  mathematical  construct 
for  determining  the  statistical  distribution  of  the  success  probability  for  a  random  system  whether 
it  be  physical  or  biological.  However,  this  distribution  alone  lacks  the  ability  to  connect  the 
system’s  response  to  an  explanatory  variable.  This  deficiency  is  remedied  by  replacing  p  with  a 
probability  calculated  from  the  cumulative  distribution  function  developed  from  the  fitting 
procedure  described  in  Section  2.3.  This  function  depends  upon  the  distribution  produced  by  fitting 
to  the  test  data.  Therefore,  it  isomorphically  maps  to  the  explanatory  variable  through  the 
distribution  function.  This  result  allows  confidence  intervals  to  be  defined  on  the  explanatory 
variable  axis.  By  doing  so,  success  probabilities  can  be  calculated  in  terms  of  the  explanatory 
variable  for  the  standard  normal  distribution  via  (2.2.2)  and  (2.2.3).  To  illustrate  this  concept,  let 
us  revisit  the  example  scenario  described  in  Sections  2.7  and  2.7.1. 

To  recap,  the  issue  under  consideration  is  determining  confidence  intervals  for  drug 
efficacy.  It  is  important  to  determine  the  ideal  dosage  required  for  an  antibiotic  to  achieve  a 
curative  blood  concentration  with  a  chosen  level  of  reliability.  The  dosage  (perhaps  expressed  in 
milligrams)  serves  as  the  explanatory  variable.  In  the  determination  of  dosage-mortality  curves, 
the  natural  logarithm  of  the  dosage  is  sometimes  used.  Sensitivity  tests  are  conducted  at  a  set  of 
discrete  dosage  levels  in  order  to  calibrate  a  normal  distribution  function  for  the  dmg’s  success.  A 
fairly  common  practice  is  to  estimate  the  drug  dosage  that  will  cure  99%  of  the  test  subjects. 
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Accordingly,  dosage  confidence  intervals  may  be  computed  for  a  confidence  coefficient  of  95%. 
Starting  with  the  theory  presented  in  the  preceding  section,  we  can  construct  this  interval. 

Upon  revisiting  equations  (2.7.1.12)  and  (2.7.1.13)  for  Model  1  or  (2.7.1.16)  and  (2.7.1.17) 
for  Model  2,  note  the  presence  of  the  term  p  ,  an  estimate  of  the  exact  success  probability.  Based 
upon  our  fit  to  the  normal  distribution,  endpoint  probabilities  plow  and  pmgh  can  be  mapped  onto 


corresponding  arguments  t Low  and  tHigh  for  the  standard  normal  distribution  ;v(o,l),  i.e. 

1 


Plow  ~ 


a/2  7T 


Lmv 

1  exp 

j 

—00 

{  2, 

dt 


(2.7.2. 1) 
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The  arguments  t Low  and  tHigh  are  functions  of  the  endpoint  probabilities  and  can  be  extracted  from 
(2.7.2. 1)  and  (2. 7. 2. 2)  via  numerical  solution  procedures.  The  associated  probits  are  related  to  the 
dosage  interval  endpoints  XLow  and  xHi  h  through  the  equations 


,  _Xlow~P.  t 
*  Low  —  ’  H‘Sh 
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XH,gh~P 
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C 2.72.3 ) 


where  p  and  a  are  the  maximum  likelihood  estimates  for  the  mean  and  standard  deviation  obtained 
by  the  techniques  described  in  Section  2.4.  On  the  other  hand,  if  a  linear  fitting  form  is  used  as  in 
Section  2.3,  the  following  forms  are  applied  in  lieu  of  (2.1 .22). 


t Low  a  ^  XLow  '■>  t High  -  a+b  XHigh 


(2.1.2A) 


The  values  tLow  and  tHigh  are  implicit  functions  of  their  corresponding  probabilities.  In  fact,  within 
the  context  of  archived  literature,  this  function  is  denoted  as  the  “Probit”.  Therefore,  (2.7.2. 1)  can 
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be  expressed  as 

tLow=  Probit  (PloJ  (2.72.5) 

A  similar  formula  may  easily  be  written  for  tmgh. 


2.7.3  Confidence  Analysis  for  the  Mean  and  Standard  Deviation 

The  content  of  this  section  constitutes  an  upgrade  to  the  corresponding  sections  in  the 
predecessor  of  this  report. [16]  The  analysis  shown  here  is  a  modified  version  of  the  derivations 
presented  in  the  DiDonato’s  excellent  report. [17]  As  is  true  for  the  maximum  likelihood  analyses 
described  above,  the  crux  of  the  probit  fit  resides  in  the  interplay  between  the  binomial  distribution 
(a  discrete  probability  model)  and  the  normal  distribution  (a  continuous  probability  model).  For 
this  reason,  when  expected  values  are  computed  for  the  log-likelihood  function,  expectation 
operators  must  be  applied  for  both  the  discrete  and  continuous  distributions.  To  assess  the  level  of 
confidence  in  the  mean  and  standard  deviation,  begin  with  the  binomial  distribution,  i.e., 


fi(s;ju,o-)  = 


^  n  ^ 


(«,— s)  S 
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(2.7.3. 1) 


where  n  =  0,l,...,n.  Recall  that  probability  p(x\ju.a)  is  given  by  the  normal  distribution  with  q  = 
1  -  /?;  note  that  x  is  the  explanatory  variable  for  the  normal  probabilities.  The  explanatory  variable 
is  not  an  active  participant  in  this  analysis.  Define  the  functions  G  and  H  as  follows. 
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G,=— [in  ./&;//,  <r)]  (2.7. 3. 2) 

O/l 

Hi  =^-[\nfi(s;Ju,cr)]  (2.7.33) 
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At  this  point,  it  is  desirable  to  take  the  expected  values  of  G  and  H.  Begin  by  treating  /  as  a 
continuous  function  of  s;  although  s  takes  on  discrete  values,  we  can  regard  s  in  the  sense  of  a 
distribution.  Hence, 


£(G,(s;//,o-)) 


\ 

Gt  (5;  (5;  ju,<r)ds 


(2.13.4) 


where  Ed  represents  the  discrete  expectation.  Substituting  for  G, 


(2.7. 3. 5) 


The  natural  logarithm  has  special  properties  under  differentiation,  i.e., 

E(G,  (s-  n,o))  =  Ed{  J  — ^ - t  f.  (« y ;  M,a)ds 
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Clearly,  we  have  that 
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(2.13.1) 


Consider  the  integrand  and  perform  the  differentiation. 
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By  differentiating  the  probabilities,  noting  q  =  1  -  p  and  collecting  terms,  we  obtain 


dju 


Gi  ^ 


v  5  y 


v^Ay 

9,_1  [fa  -s)q-sp]^~ 

OjU 

(2.7.3.10) 

qEAsh>  [Ed (n  —  s)q  —  Ed (s)p]^~ 

O/l 

(2.7.3.11) 
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By  nature  of  the  binomial  distribution,  both  (n  -  s)  and  S  are  discrete  random  variables.  Moreover, 

Ed ini  ~s)  =  nP  ;  Ed  W  =  ni 1 1  (2.7.3. 12) 

Substitution  in  (2.7.3.10)  yields 

Ed[^-\  =  (n‘  1  P",P^  cin'9~'  [n,Pcl  ~  nicip]^~  =  0  (2.7.3.13) 

yd/u)  y  s  J  dju 

By  backing  these  results  into  (2.7. 3.7),  we  can  conclude  that 

E(G,.)  =  0  (2.7.3.14) 

It  is  also  necessary  to  determine  the  expected  value  of  Hh  this  calculation  is  easily  accomplished 

because  the  only  difference  is  that  the  derivative  of  fi  is  taken  with  respect  to  <J  instead  of  // . 

The  process  is  the  same,  and  the  result  corresponding  to  (2.7.3.14)  is 

Ed[^L  1  =  Pn,P~l  q^fopq  ~  niqp\~~  =  0  (2.7.3.15) 

'ydcr  J  ys  J  ocj 

Hence, 

EiH^  0  (2.7.3.16) 

Mirroring  the  maximum  likelihood  analysis,  we  the  variances  for  Gt  and  H , ,  then  the 
covariance  of  Gt  and  H , .  Begin  with  the  variance  of  Gh  i.e., 

Var(G()  =  <j2Gj  (2.7.3.17) 

By  definition  with  the  use  of  (2.7.3.14), 

<  =  £(g(2)-(£(G,))!  =  fi(G,2)  (2.7.3.18) 

dc/b  JG(2(^;^.o-)/,(5;//,(j)</5  (2.7.3.19) 

s 

Substituting  for  G  and  then  evaluating  the  derivative  of  the  logarithm. 
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^(g,2  )  =  J  ^ln  fi  (5;  fi,  cr)  f  (5 ;  H,  o)ds 


(2.7.3.20) 


.(^2)  f  1  dfi(s\ju,c j)  (  ^ 

HG'  i=  777 - \  -  fi[s-^(7)ds 

s  fi  {s',M,°)K 


(2.1.3.21) 


Thus,  the  initial  result  is 


7)=f- 


i  &U 


f\dM 


(2.7.3.22) 


Now  consider  the  partial  derivative  with  respect  to  /J.  of  E  (G, ) ,  i.e., 


cEiG,)  _  d 
dju  dju 


J  Gi(s‘,ju',cr)fi(s',ju,<j)ds 


(2.7.3.23) 


Recall  (2.7.3.14)  and  interchange  the  order  of  differentiation  and  integration;  note  that  the  domain 
of  integration  is  finite,  and  the  differentiation  is  with  respect  to  a  parameter  not  the  variable  of 
integration. 


0  =  [  —  [G,  (5;  ju,  cr)fI  (5;  ju,  a)]ds 
Jd/u 


(2.7.3.24) 


By  applying  the  product  rule, 


r  dGt(s',  ,  x  (  \dft(s; n,<r) 

0  =  I - - - -fi  (5;  //,  a)  +  G,  (5;  //,  a) - - - ds 

J  dju  dju 


(2.7.3.25) 


Using  (2.7. 3. 2)  to  substitute  for  G ,, 


j  d^2  )  j  ^ u  dju 


(2.7.3.26) 


Rearrangement  along  with  the  properties  of  the  logarithm  yields 


1  f  df.(s;/i,er) 


f{  fy 


ds  =  -  j  -f  {s',  ju,a)ds 


(2.7.3.27) 
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Upon  substitution  into  (2.7.3.22)  we  learn  that 


E(Gi2 )  =  -J  dJ^fL  ft  (5;  P'  °)ds 


(2.7.3.28) 


By  definition  of  the  expectation, 


e(g?)=-e 


rd2\nf^ 

d/u 2 


(2.7.3.29) 


The  variance  of  Hi  is  computed  via  a  similar  analysis;  the  only  real  change  is  that  we 
differentiate  with  respect  to  a  instead  of  [i.  Therefore, 


al=E(Hf)=-E 


rd2\nf^ 
v  5(j2  j 


(2.7.3.30) 


It  is  also  necessary  to  determine  a  form  for  the  covariance  of  G,  and  H,.  This  expression  is 
derived  as  follows. 


Cov(G„  H,)  =  <r  =  jain/  ain  f,  a)ds 

-  dju  da 


(2.7.3.31) 


As  before,  the  properties  of  the  natural  logarithm  are  applied. 


W=J 


1  a/,  1  a/, 


An  =  —  - - - 


fi  5/v  fi  da 


fi(s;Ju,a)ds  =  J 


1  dfi  dfi 

ft(s\n,a)  d/u  da 


ds 


(2.7.3.32) 


Recall  (2. 7. 3. 2)  and  differentiate  it  with  respect  to  cr,  i.e., 

E{Gt )  =  |  G,/;  (s;  //,  a)ds 

S 

Using  (2.7.3.14)  and  interchanging  the  order  of  integration  and  differentiation, 

o  =  MS)  =  A  J  Gifi(s\  a)ds  =  f  ^  *  -)) 
da  da J  J 


(2.7.3.33) 


da 


-ds 


(2.7.3.34) 


By  applying  the  product  rule, 
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Substituting  for  G,,  we  obtain 


0  =  {  9  ft  (, s ;  //,  cr) +  J 


d/uda 


8 In/,  5/, 

•J  0//  0CT 


Again,  the  logarithm’s  properties  are  applied  to  obtain 

In  ft  ,/ _ ,  f  1  5/;  0/. 


0  =  J  h  fi(s;/J,cr)ds  +  | 

J  OjUOCT  J 


Hence, 


1 


1  5/,  5/, 


ft  dfu  da 


ds 


r  0  In  /.  /  n 

a  j1  =  -  - -  fAs ;  ju,a)ds 

fi  dju  da  •;  djuda 


Substituting  this  result  into  (2.7.3.32),  we  conclude  that 

02  In  /, 


=  -  j  /,  (^;  //,  <y)ds  =  -E 

J  djuda 


fd2\nf^ 

djuda 


(2.7.3.36) 


(2.7.3.37) 


(2.7.3.38) 


(2.7.3.39) 


having  used  the  definition  of  the  expectation.  Equation  (2.7.3.32)  demonstrates  that  this  covariance 
is  completely  symmetric,  i.e., 


Cov  (G, ,H.)  =  aGH  =  (JHjGi  =  Cov  {Hi ,  Gi )  (2.7.3.40) 

This  development  is  quite  lengthy,  but  at  this  point,  we  can  begin  to  roll  up  the  results 
shown  above  into  a  form  pointed  toward  a  confidence  assessment.  Recall  (2.7.3. 1),  the  factor 
of  the  likelihood  function;  it  follows  that  the  entire  likelihood  function  is  written  as 

N 

F  =  Hfi{s-,ju,a )  (2.7.3.41) 

(=i 

Recall  that  the  index  i  ranges  over  the  A  levels  taken  on  by  the  explanatory  variable.  Accordingly, 
the  log-likelihood  function  is  expressed  as 
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(2.7.3.42) 


L  =  lnF  =Yjlnfi(s;Ju,o) 


1=1 


Equation  (2.7.3.42)  is  an  important  result  because,  via  the  properties  of  the  logarithm,  the  product 
of  probabilities  in  the  likelihood  function  has  been  transformed  into  a  sum  of  N  log-likelihood 
factors.  Using  (2. 7. 3. 2)  and  (2. 7. 3. 3),  partial  derivatives  (taken  with  respect  to  the  fitting 
parameters)  of  the  log-likelihood  function  can  be  written  as 

dL  _  A  d  In  /( (.v;  p,  o) 


y 

dju  !=1  d/i 


TdGi(s\ti,a) 


i= 1 


dL  ^dlnf^/1,0)  ^  u  (  ^ 


1=1 


These  derivatives  can  be  organized  in  the  form  of  a  vector 


N  \ 


'  dL  dL  ] 

f 

y  d/j  do  j 

V 

(2.7.3.43) 


(2.7.3.44) 


(2.7.3.45) 


For  large  sample  sizes,  that  being  defined  by  the  used  at  each  level  of  the  explanatory  variable, 
the  vector  (2.7.3.45)  is  distributed  in  multinormal  distribution  TV (o,  A  1  ]  where 


A-1  = 


E\ 

e\ 


'  d2L ^ 

J7 

'  d2L  ' 

l  dF2  J 

L/ 

dju  do  J 

'  d2L  'l 

'  d2L ' 

v  dju  do  j 

l  3(72J 

Z' 

1=1 

N 

Z' 

i=l 


Gi  GjHj 

i=l  i=l 

N  N 

_2 


GjHj 


z< 

i=l 

N 

Z' 


H, 


i= 1 


V  A*, 


(2.7.3.46) 


A-1  is  denoted  as  the  covariance  matrix.  [17]  Using  this  information,  the  departure  of  the 
parameter  estimates  from  their  true  values  can  be  expressed  in  the  form  of  a  Taylor  cast  in  two 
variables,  i.e., 

=  YLU.°'o)+TTU.ffi)(/'-//o)  +  77L(/'i.ffi)(o;-ffo)  (2.7.3.47) 

0/1  O/l  O/l  O/lOO 
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^L{p,cr)  =  ^L  (a,  ,  cr0 )  +  (/T ,  cr2 ) (p  -  Mo )  +  |~t  (/T  >  <*i ) (&  -  <r0  )  (2.7.3.48) 

0(7  0(7  OjUOa  0(7 

These  expressions  are  presented  as  equalities  in  deference  to  the  standard  form  for  the  Taylor 
series.  Instead  of  expanding  the  second  order  partial  derivatives  around  the  point  (//0,<70)?  they 
are  centered  at  two  points  ( // \ ,  cr, )  and  (p2 ,  <r2 )  that  lie  along  the  line  segment  between  (p,  a)  and 

(/./„ ,  <T0  ) .  Although  these  points  are  not  known  a  priori,  they  both  converge  to  (p.  a)  with 
probability  one.  Moreover,  the  elements  of  (2.7.3.46)  converge  to  the  expected  values  with 
probability  one. [17]  At  (p.a),  the  first  order  partial  derivatives  vanish,  so  (2.7.3.47)  and 


(2.7.3.48)  are  rewritten  as 


~p~  (Ao  ’  ) : 

op 


8  L  (A, ,  o', )(//  -  Ao ) - (A, , )(^  - )  (2.7.3.49) 


dp 


dpdcr 


dL  d2L  d2L 

T—  (Ao ,<r0 )  =  - (a2 >  ^2 ){p  -  Mo ) ■ -  (a2 ,o-2){&-(7 o )  (2.7.3.50) 

da  dp  da  da 

Finally,  at  this  point  of  the  development,  we  can  address  the  level  of  confidence  in  the 

estimates  /7  and  cr  .  To  do  so,  it  is  instructive  to  write  the  linear  system  (2.7.3.49)  and  (2.7.3.50) 

in  matrix-vector  form.  Define  vectors  as  shown  below: 


9  = 


M  ~  Mo 
a-a0_ 


/L  = 


dp 

(/4)  ’  G"o  ) 

. da 


(2.7.3.51) 


The  (2.7.3.49)  and  (2.7.3.50)  can  be  written  as 

A  =  Al*0 


(2.7.3.52) 


or 


6  =  A  •  X 


(2.7.3.53) 
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The  vector  0  has  an  approximate  bivariate  multinormal  distribution,  that  is  n(0,a).  The  matrix 
A  is  denoted  as  the  asymptotic  variance/covariance  matrix.  To  fix  a  level  of  confidence  to  the 
estimate  (jj,  cf) ,  consider  the  quadratic  form 

Q  =  6T*A~x*e  (2.7.3.54) 

For  a  vector  X  with  distribution  n(v,v)  ;  the  associated  joint  density  function  f(x)  has  the  general 
form 

/(*)  =  (^pfqexp[“i(X  ~v'l'V *(X  -v)r]  (2.7.3.55) 

2 

This  density  function  has  a  J"  (Chi-squared)  distribution  with  m  degrees  of  freedom  where  m  is 
the  number  of  components  in  the  random  vector  X.  [14]  For  our  specific  concern  (2.7.3.54),  there 
are  two  degrees  of  freedom:  for  the  Golub-Grubbs  model  and  (a.  /?)  for  the  Garwood 

2 

model.  For  the  former  model  with  Q  provided  by  the  %  "  distribution,  (2.7.3.54)  can  be  expanded 
to  obtain 

(m  -  mY  +  2 V (m  -  z7)^  -  °)  +  A«o (o-  -  °Y  =  zlr  (2.7.3.56) 

In  this  equation,  the  calculation  is  centered  at  the  parameter  pair  estimate  (jj,  rf) ,  and  so  are  the 

matrix  elements  A  ,  A^a  and  Aoa .  The  locus  defining  the  boundary  for  the  confidence 

assessment  is  given  by  the  ordered  pairs  (ju,a),  a  continuum.  The  position  of  this  locus  is 

2 

controlled  by  the  term  on  the  right  hand  side,  the  tail  of  the  /  ’  distribution  for  confidence  level 

2 

y.  The  /  term  is  given  by 


(2.7.3.57) 


Distribution  A 
45 


Hence, 


Xly 


-21n(l-Pr(j2  <1-^)) 


(2.7.3.58) 


As  an  illustration,  for  a  confidence  level  of  y  =  0.95  ,  1  -  y  =  0.05  ,  the  value  of  X\~r  is  -0.1025.  If 

one  cares  to  apply  some  analytic  geometry,  (2.7.3.56)  can  be  shown  to  be  an  ellipse  in  the  plane. 
The  ellipse  can  be  rotated  into  a  standard  position,  and  the  lengths  of  its  semi-major  and  semi¬ 
minor  axes  can  be  calculated.  Doing  so  is  all  well  and  good,  and  it  is  even  instructive,  but  it  adds 
complication  through  the  rotation.  Instead,  we  can  simply  solve  the  quadratic  equation  for  one 
parameter  of  the  other.  If  we  solve  for  G  ,  we  obtain 


a  =  a  - 


Xi- 


-  + 


f  A2 

^ H<T 


w 


(a-a)2 


(2.7.3.59) 


y  GG  OG  J 


To  plot  this  ellipse,  we  must  determine  the  domain  of  the  independent  variable  /u  .  To  do  so,  the 
argument  of  the  radical  must  be  non-negative,  i.e., 


xL 


+ 


(  A2  A  ^ 

A 2  A 

\  ^GG  **~GG  J 


(//-//  )2  >0 


(2.7.3.60) 


The  endpoints  of  the  interval  for  /J  are  determined  by  solving  (2.7.3.60)  for  the  equality.  These 
endpoints  are: 


An 


A 


zLA 


y  go 


Ka  -  KAo° 


(2.7.3.61) 


Amax  =  A  +  • 


ZlyA 


y  gg 


(2.7.3.62) 
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Therefore,  //  e  |//mm  ,  //max  J  •  Within  this  domain,  (2.7.3.59)  defines  two  branches  of  the  ellipse  in 


CF  that  are  easily  plotted  so  long  as  the  radical  is  made  to  vanish  at  the  endpoints.  Note  that  floating 
point  arithmetic  tends  to  create  extremely  tiny  negative  numbers  under  the  radical.  The  argument 
must  be  set  to  zero  in  these  cases. 

After  the  pages  and  pages  of  mathematics  required  to  reach  this  point  in  our  development, 
it  is  of  great  importance  to  interpret  the  meaning  of  this  confidence  assessment.  The  confidence 
ellipse  exists  in  the  parameters  plane  whereas  the  confidence  interval  consists  of  only  2  discrete 
points.  DiDonato  [17]  claims  that  the  true  parameter  pair  lies  within  the  ellipse  with  a  probability 
of  7.  This  claim  lies  in  disagreement  with  the  frequentist  interpretation  of  probability.  At  the 
conclusion  of  a  single  experiment  consisting  of  n  levels  in  an  explanatory  variable,  a  single 
parameter  pair  is  calculated;  then  the  elliptical  locus  is  calculated.  These  calculations  are  entirely 
deterministic;  there  is  nothing  random  about  them.  The  only  way  that  randomness  is  restored  is 
through  repeated  experiments.  That  is  to  say,  the  entire  experiment  must  be  repeated  say,  M  times. 
Each  repeat  of  the  experiment,  becomes  a  separate  realization  of  the  experiment.  Moreover,  each 
repeat  of  the  experiment  may  use  different  levels  and  number  of  levels;  it  may  also  use  different 
numbers  of  trials  at  each  level.  In  this  sense,  randomness  is  restored  to  the  ensemble  of  experiments 
and  to  the  locus  of  points  constituting  the  confidence  ellipse.  With  this  interpretation,  if  10 
experiments  are  performed  with  the  subsequent  calculation  of  90%  confidence  ellipses,  then  9  of 
the  10  ellipses  will  contain  the  true  parameter  pair.  For  any  of  the  individual  experiments,  then  the 
probability  that  the  experiment’s  confidence  ellipse  contains  the  true  parameter  pair  is  either  one 
or  zero.  Taken  by  itself,  the  single  experiment’s  ellipse  only  provides  the  analyst  with  a  sense  of 
the  extent  of  scatter  in  the  parameter  pair.  Scheffe  [17]  echoes  this  interpretation  while  also  citing 
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the  former  claim,  but  that  latter  interpretation  is  rigorous  from  frequentist  theory  and  is  adhered  to 
in  this  report. 


2.8  Measures  of  Variance 

In  this  section,  we  address  measures  of  variance  associated  with  our  fitting  procedure.  As 
it  happens,  estimates  of  the  mean  /u  and  standard  deviation  a  for  the  probability  distribution,  have 
their  own  variances  Var(w)  and  Var(cr).  Also,  we  can  imagine  that  (i  and  a  vary  in  a  correlated 
manner  in  the  probability  distribution.  This  eventuality  implies  that  the  covariance  of  /./  and  a  may 
be  non-zero.  As  it  happens,  these  estimated  variance  values  are  given  as  entries  in  Avar ,  the 
asymptotic  variance-covariance  matrix,  the  inverse  of  the  information  matrix.  The  information 

matrix  Ainfo  is  given  by 


Anfo=-£ 


d2L  d~L 


dju 2 

d2L 


djiida 

d2L 


(2.8.1) 


djuda  da 2 

The  asymptotic  variance-covariance  matrix  is  the  inverse  of  the  information  matrix,  i.e., 


Avar  = 


-E\ 

-E 


' d2L ) 

-  F 

'  d2L  ' 

KdE2  J 

ydjuda  J 

“i-i 


'  d2L  ) 

-  F 

"a2L" 

Kdjudcr  j 

J-j 

{da2} 

(2.8.2) 


The  second  partial  derivatives  are  calculated  during  the  estimation  procedure,  as  in  Section  2.6.  It 
is  easy  to  compute  their  expected  values  by  using  properties  of  the  binomial  probability 
distribution.  The  elements  of  the  asymptotic  variance-covariance  matrix  are  interpreted  as  follows. 
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Avar  - 


Var(//)  Cov(pi,  a) 
Co  v(pi,a  )  Var(  a) 


(2.8.3) 


This  matrix  contains  the  desired  variance  estimates  for  the  parameters  produced  the  fitting 
procedure. 

On  the  other  hand,  suppose  that  the  linear  fitting  form  is  used  as  in  Section  2.3,  i.e., 

t  =  a  +  fix  (2.8.4) 

The  information  matrix  elements  cast  in  terms  of  a  and  [>  can,  by  transformation,  be  expressed  in 
terms  of  pi  and  a.  The  information  matrix  expressed  in  terms  of  a  and  [1  is  written  as 

32L  82L 


Anfo  - 


da2 

d2L 


dad/3 

d2L 


(2.8.5) 


dad/3  d/3 

Recall  that  the  transformation  (2.3.2)  has  the  general  form 

«  =  «(//,  a) ;  J3  =  j3(/j,cr) 

By  using  the  chain  rule  for  partial  differentiation,  we  observe  that 


d  _da  d  +d(3  d 


(2.8.6) 


dpi  dpi  da  dpi  d/3 


d  _  da  d  d/3  d 
da  da  da  da  d/3 


(2.8.7) 


(2.8.8) 


It  is  a  tedious  mathematical  exercise,  but  the  repeated  application  of  (2.8.7)  and  (2.8.8)  leads  to 


d2L 

dpi2 


da  d  L  d/3  d2L 


dpi  da2  dpi  dad/3 
da  d2L 

+ ' 


da  dL  d2a 

- 1 - x 

dpi  da  dpi “ 


d/3  d2L 
dpi  dad/3  dpi  d/31 


■  + 


d(3_  dL  d2 (3 
dpi  d/3  dpi 2 


(2.8.9) 
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82L 

djuda 


(2.8.10) 


da  d  L  dfi 


d2L 


da  da 2  da  dadfi 

da  d2L  dfi  d2L 


da  dL  d2a 

- 1 - 

d/j  da  djuda 


+ 


-  +  - 


da  dadfi  da  dfi~ 


dfi_  dL  d2fi 
dju  dfi  djuda 


d  L 

da2 


da  d~L  dfi 


d  L 


da  da2  da  dadfi 

da  d2L 

+ ' 


da  dL  d2a 

- 1 - 7 

da  da  da " 


dfi  d2L 
da  dadfi  da  dfi2 


■  + 


dfi_  dL  d2fi 
da  dfi  da2 


(2.8.11) 


Equations  (2.8.9)  through  (2.8.11)  are  the  elements  of  Awo,  and  they  are  easily  computed  with  the 
use  of  the  following  coordinate  transformation  derivatives. 


da 

dju 

1 

9 

G 

da  ju 

da  a2  ’ 

fie  l-gf 

T:  |  os 
II 

o 

dfi  _  1 

da  a2 

(2.8.12) 

52a_0 

d2a 

1 

d2a 

=  _2jJ_ 

(2.8.13) 

d^2 

djuda 

a2 

da2 

a 3 

d2fi  _ 
dju 2 

0;  ^ 
djuda 

(O  CO 

o 

II 

2 

—  F 

a 

(2.8.14) 

Equations  (2.8.9)  through  (2.8.14)  provide  the  required  information  matrix  entries  needed  in  order 
to  compute  the  asymptotic  variance-covariance  matrix.  For  purposes  of  calculation,  (2.8.9)  through 
(2.8.11)  can  be  simplified  by  noting  that  dL/ da  and  dL/ dfi  are  zero  at  the  point  of  maximum 
likelihood. 


2.9  Logit  Fitting  Procedure 

The  probit  analysis  discussed  earlier  in  this  report  is  based  upon  the  normal  distribution.  In 
some  cases,  it  is  desirable  to  use  a  distribution  that  performs  a  bit  differently  in  the  tail  region.  The 
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logistic  distribution  is  often  applied  for  linear  fitting.  This  distribution  is  computed  as  the  natural 
logarithm  of  the  odds  ratio,  i.e., 


logit(p)  =  in 


1  ~P 


=  a  +  fix 


(2.9.1) 


In  this  formula,  p  is  the  success  probability,  so  the  odds  ratio  is  clearly  present.  The  explanatory 
variable  is  x  while  a  and  fi  are  the  unknown  parameters.  Rearrangement  of  (2.9.1)  allows 


calculation  of  the  probability;  specifically, 


P{x ;  a,  fi) 


exp  (a  +  fix) 
l  +  exp(a  +  fix) 


(2.9.2) 


Note,  unlike  the  probit  method,  the  computation  of  this  probability  does  not  require  the  use  of 
numerical  quadrature.  The  probability  density  function  is  calculated  by  differentiating  (2.9.2)  with 
respect  to  x.  Hence, 

dp{x\a,fi)  =  pzvp{a  + fix) 
dt  [l  + exp  («  +  /?*)]" 

To  elucidate  some  of  the  properties  of  this  distribution,  we  can  substitute  t  =  -a  -  fix  then  rewrite 

(2.9.2)  as 


f 


It 


P 

~Pj 


=  -t 


(2.9.4) 


exp(-f) 

1  +  exp(- 1) 


(2.9.5) 


The  complementary  probability  is 


1 

1  +  exp(-f) 


(2.9.6) 
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For  application  of  this  distribution,  t  can  be  substituted  for  ~t .  However,  for  the  proof  below, 
(2.9.5)  must  be  left  unaltered.  The  expected  value  for  this  distribution  is  zero,  i.e.,  e(i)  =  0  .  This 


property  is  shown  by  following  [18].  By  definition, 


t exp(- 1 ) 


Eif)  =  J  =  J  y - 

L  -oo  [l  +  exp(- 1 


)1 


-dt 


This  integral  can  be  rewritten  as  follows. 

t exp(- 1 ) 


£(<W 


oo 

dt  +  \~\ 


t exp(- 1 ) 


I  [l  +  exp(-  0] 2  o  [l  +  exp(—  t )] 


-dt 


In  the  leftmost  integral,  make  the  substitution  t  =  —y  . 


£W=J 


y  ex 


P(v) 


dy  +  \} 


texp(-t) 


o  [l  +  exP(.v)]2  '  o  [•  +  exp(-  /)] 


-dt 


By  expanding  the  quadratic  term,  we  obtain 

-yexp(y) 


£W=JT 


p  +  2exp(v)  +  exp(2v)  '  ]  1  +  2exp(-^)  +  exp(- 2l) 


t exp(- t ) 


-dt 


(2.9.7) 


(2.9.8) 


(2.9.9) 


(2.9.10) 


If  we  multiply  the  integrand  of  the  left  integral  by  exp(-  v)/exp(-  y)  and  the  integrand  of  the 
right  integral  by  ex  p(0/  exp(r) ,  we  have  that 


£(/)=  f - 1 - dz+  f - - - dz  (2.9.11) 

J0  exp(-z)  +  2  +  exp(z)  J0  exp(z)  +  2  +  exp(-z) 

Since  z  is  a  dummy  variable  of  integration,  these  two  integrals  are  of  the  same  magnitude  but 
opposite  sign;  thus, 


E(t)  =  0 

For  the  linear  form  shown  in  (2.9.1),  (2.9.11)  implies  that 


E(a+Px)  =  ~ 


(2.9.12) 


(2.9.13) 
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The  variance  for  this  distribution  can  also  be  calculated.  By  definition,  and  with  the  use  of  (2.9. 12), 

Var  (t)=E(t2}-(E(t))2=E(t2)  (2.9.14) 

This  expectation  may  be  evaluated  in  the  same  manner.  Observe  that 

F(f2)—  f  (2  exP(~  0  dt-  f  rexp(-t)  ,,  ,  f  t2  exp(-t) 


[l  +  exp(-  /)] 2  i  [l  +  exp(- f )] 2  J0  [l  +  exp(-  /)] 


dt  +j*  -j 


-dt 


(2.9.15) 


Again,  the  substitution  t  =  -y  is  employed  in  the  integral  over  negative  I,  i.e., 


q  .2 


r  t°exp(-t)  r  (- y)  exp(y)  r 

r -  ;  L  dt  =  — T^idy  = 

,[l  +  exp(-t)]  ^  h  j. J 


y  ex 


p(y) 


Simplification  shows  that 


[l  +  exp(y)]2 


[exp(y)(exp(-y)  +  l)]: 


-cly 


(2.9.16) 


0  .2 


1 


t  ex 


p(-0 


-I,  [l  +  exp(-  /)] 


y~  ex 


p(-  y) 


[(exp(-.v)+l)] 


-dy 


(2.9.17) 


This  integral  is  the  same  as  the  remaining  integral  over  positive  t  in  (2.9.15),  so 

t2 exp(- 1 ) 


49=2/ 


o  [l  +  exp(- /)] 


-dt 


(2.9.18) 


To  evaluate  this  integral,  let  y  =  exp(-r),  with  this  substitution  the  integrand  in  (2.9.18)  becomes 


.V 


(i+A 


Z(-r«/ 


(2.9.19) 


The  series  expansion  appearing  in  (2.9.18)  is  derived  by  applying  polynomial  division  to 
1  /(I  +  y)2  and  then  multiplying  by  the  z  remaining  in  the  numerator.  Therefore, 

<P(-  <)  _ 


ext 


[l  +  exp(-f)]2 


Z  (-  •)"' n  (exp(-  *))"  =  Z  (-  O"1  n  exp(~ nt) 


(2.9.20) 


Substituting  this  result  into  (2.9.17), 


/  \  00  °° 

E(t2 )  =  2^  (- 1)"-1  nj  t2  exp(-  nt)dt 

n=  1  0 


(2.9.21) 
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Integration  by  parts  can  be  used  to  evaluate  the  integral  in  (2.9.20),  i.e., 


Jr  exp(-  nt)dt 


With  this  result,  we  learn  that 


4!)=2Z(-ir«f4l=4Z^ 


n= 1 


\n  J  „=i  n 


The  convergent  series  in  (2.9.22)  evaluates  to  a  value  of  n2 1 12,  so 

44 


n 

T 


If  we  return  this  result  to  the  original  logistic  form  (2.9.3),  we  see  that 

Var(a  +  fix)  =  E(a  +  fixf 


fi2n 


Equations  (2.9.12)  and  (2.9.24)  convey  basic  properties  of  the  logistic  distribution, 
employed  for  analyzing  sensitivity  test  data.  The  likelihood  equation,  applied  again,  is 


'  =  n 


n  rn  \ 


S:  n.:  —S: 


Pi1  <7, 


so  the  log-likelihood  function  is 


N  (n\ 


\SiJ 


+  iV,  In  Pi  +  {nj  -  si  )ln  qt 


The  log-likelihood  function  is  maximized  by  computing  the  derivatives 

dL  4  dP,  !  (n,-s,  )dch 


da  pt  da  qi  da 
8L  ^  dPi  |  (n^sfidcji 


dfi  i= i  Pi  dfi  q,  dfi 


For  this  application,  the  success  probability  has  the  form 


(2.9.22) 

(2.9.23) 

(2.9.24) 

(2.9.25) 
It  can  be 

(2.9.26) 

(2.9.27) 

(2.9.28) 

(2.9.29) 
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(2.9.30) 


Pt  = 


exp  (f,.) 


l  +  exp(ff) 

For  readers  who  have  read  through  previous  discussions  in  this  report,  the  terms  in  the  two 
equations  above  will  appear  very  familiar,  and  indeed,  they  are.  For  ti=a  +  fixi ,  we  can  use  the 
derivatives  provided  in  Section  3.  It  is  particularly  important  to  note  that  dp;  I  dti  =  zt ,  i.e., 

exp(0 


(l  +  exp(r . )) 2 

With  these  derivatives  in  mind,  (2.9.27)  and  (2.9.28)  are  rewritten  as 

dL 


da 


Z*,-— -(ni~si)— 


i= 1 


dL  Zj  /  \  Zj 

-zxy  =  Ls,x,  —  \ni  ~  si )xi  — 

dp  i= i  Pt  qt 


The  attendant  second  partial  derivatives  are  expressed  as 

d2L 


dec  /=,  {Pi  J  yep  J 


d~L 


fz  A 


~  Xn=HSiXi  —  ~(ni~Si)Xi 
dadp  tT  ^  p( ) 


V 

\d,  j 


(2.9.31) 


(2.9.32) 


(2.9.33) 


(2.9.34) 


(2.9.35) 


d~L 


dp 


f  x 


2  -  ^SiXi 


\PiJ 


-  (ni  -  si)xf 


f  X 

Ll 

U, 


(2.9.36) 


The  prime  '  notation  represents  the  ordinary  derivative  d  I  dt{ .  The  derivatives  of  ratios  occurring 
above  are 


(  X 
Zi 

_  d 

(z^ 

UJ 

dtt 

l  Pu 

PiZi-PiZ 


i 


Pi  Pi 


(2.9.37) 
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(2.9.38) 


f  \ 

_  cl 

(  \ 
h. 

■Ml  J 

dti 

■Ml  y 

ch  Zj  ~  ch  Zj 
2 

Qi 


Equations  (2.9.32)  through  (2.9.38)  may  be  used  in  a  Newton  iterative  scheme  as  described  in 
Section  2.5. 

The  logit  procedure  (as  described  above  and  upon  convergence)  renders  an  estimate  of  the 
parameters  a  and  ft .  Naturally,  questions  arise  as  to  the  certainty  of  these  estimates.  Given  that 
the  framework  for  the  success  probability  still  resides  in  the  binomial  distribution,  the  concept  of 
the  confidence  interval  for  the  success  probability  is  retained  here  altered  only  by  the  application 
of  the  logistic  probability  formulas  provided  above. [4]  As  in  the  case  of  the  probit  model  (based 
on  the  normal  distribution),  it  is  desirable  to  obtain  a  bivariate  confidence  relationship  between  a 
and  (3 .  Given  that  the  logit  method  utilizes  the  same  mathematical  underpinnings,  i.e.,  the 
information  and  asymptotic  variance/covariance  matrices,  one  might  suspect  that  the  same 
methodology  leading  to  the  computation  of  the  confidence  ellipse  can  be  applied.  It  is  a  promising 
idea  indeed,  but  a  troublesome  question  arises.  The  theoretical  development  described  in  Section 
2.7.3  is  based  not  only  on  the  normal  distribution  but,  more  importantly,  on  its  connection  to  the 
chi-squared  distribution  that  governs  the  uncertainty  in  parameters  //  and  a  via  a  special 
quadratic  form. [14]  This  quadratic  form  involves  both  the  information  matrix  and  the  parameter 
increments.  The  logistic  distribution  differs  from  the  normal  distribution,  so  the  connection  to  the 
chi-squared  distribution  is  not  immediately  evident.  Fortunately,  research  has  shown  that  the 
logistic  distribution  correlates  well  enough  with  the  normal  distribution  to  assert  that  the 
uncertainty  for  the  bivariate  parameters  is  distributed  chi-squared. [19,20].  For  this  reason, 
confidence  ellipses  can  be  calculated  for  logit  parameter  estimates  in  the  same  way  as  for  probit 
estimates. 
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3  TEST  PROBLEMS  AND  RESULTS 


In  the  preceding  sections,  some  of  the  probability  theory  supporting  the  analysis  of 
sensitivity  test  data  has  been  described.  The  mathematics  employed  is  quite  complicated  and  much 
time  is  required  in  order  to  gain  a  working  level  of  knowledge  in  this  discipline.  Mired  in  theory, 
it  is  easy  to  lose  touch  with  the  practical  aspects  of  the  calculations.  For  this  reason,  this  chapter 
of  the  report  presents  the  setup  and  documents  the  results  associated  with  a  series  of  basic  test 
problems.  To  enhance  understanding  and  permit  brevity,  we  have  selected  problems  from 
Garwood  [11],  Dixon/Mood  [6]  and  Collett  [4].  These  problems  are  useful  in  testing  the  algorithms 
discussed  above  in  Chapter  2.  In  and  of  themselves,  these  problem  sets  constitute  a  reasonable 
practicum  for  didactic  purposes,  but  they  do  not  address  all  of  the  potential  scenarios  for  sensitivity 
testing.  That  noble  goal  remains  on  the  horizon,  but  perhaps  it  can  now  be  glimpsed  in  clearer 
focus. 

3.1  Test  Problem  1:  Antibiotic  Efficacy 

This  example  addresses  a  series  of  serum  testing  trials  associated  with  an  anti-pneumonia 
drug.  Selected  doses  of  the  drug  are  administered  to  groups  of  mice.  The  data  is  presented  as  a 
series  of  five  binomial  trials,  but  the  data  can  also  be  represented  in  terms  of  a  sequence  of  binary 
trials  (each  treated  mouse  constitutes  a  binary  trial  at  a  particular  dosage  level).  The  data  is 
provided  in  Table  1,  and  the  drug  dosage  is  measured  in  cubic  centimeters. [11]  Each  individual 
dosage  Xd  is  administered  to  40  mice,  and  after  being  infected  with  pneumonia,  the  number  of 


Table  3  Anti-Pneumonococcus  Serum  Test  Data 


Dosage  (Xd  ,  cc) 

Successes  out  of  40 

0.000625 

-2 

33 

0.00125 

-1 

22 

0.0025 

0 

8 

0.005 

1 

5 

0.01 

2 

2 
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Figure  2  Normal  Curve  Linear  Fit  Form  for  Test  Problem  1 
expired  mice  is  counted  and  recorded  (an  expired  mouse  is  denoted  as  a  success  in  this  context). 

In  this  case,  it  is  advantageous  to  transform  the  dosage  levels  via  the  natural  logarithm.  The  dosage 

levels  shown  in  Table  3  are  cleverly  selected  to  promote  logarithmic  mapping.  Let  XD  =  In  (X D ) ; 
then  the  transformed  dosages  may  be  written  as 

X”D=(n- l)ln(2)  +  ln(Xi),  n  =  2,...,5  (3.1.1) 

Garwood  [11]  applies  a  final  transformation  to  the  data  in  order  to  center  the  distribution  at 
XD  =  0.  If  there  are  an  odd  number  n  of  data  (as  in  this  case),  then  let  m  be  the  natural  number 
with  n  =  2 in  + 1 .  The  distribution  is  centered  on  the  discrete  interval  [-  in.  m\  with  the 
transformation 


vm  — 

A  D  — 


f  ~Y  m  yl  ^ 

A  D  A  D 

*  c  ^  1 

K-Kj 


-i 


m . 


m  =  1  ,...,5 


(3.1.2) 
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Figure  3  Probability  Distributions  for  Test  Problem  1 

To  perform  the  fitting  procedure,  Garwood’s  implementation  of  the  method  of  maximum 
likelihood  (Section  2.3)  is  employed  with  the  linear  form  a  +  fix  and  the  normal  distribution.  To 
start  the  iterative  procedure,  the  initial  estimates  for  a  and  fi  are 

a  =  2;  p  =  \  (3.1.3) 

The  iterative  scheme  converges  quickly  to  the  final  values  of 

a  =  0.554428;  /?  =  0.676081  (3.1.4) 

These  values  agree  extremely  well  with  the  archived  estimates  in  Garwood.  The  linear  fit  for  this 
test  problem  is  shown  in  Figure  2,  and  the  probability  density  and  cumulative  distribution  functions 
are  shown  in  Figure  3.  In  accordance  with  the  binomial  model  (Model  1),  the  95%  confidence 
interval  endpoints  in  terms  of  the  logarithmic  dosage  for  probability  estimates  are  provided  in 
Table  4.  Confidence  interval  (95%)  endpoints  computed  using  Model  2,  the  standard  error  method, 
are  provided  in  Table  5.  Also,  the  95%  confidence  ellipse  for  the  parameter  estimates  for  a  and 


Distribution  A 
59 


Table  4  95%  Confidence  Interval  Endpoints  (Model  1)  for  Test  Problem  1 


Left  Endpoint 

In 

Probability 

Estimate 

Right  Endpoint 

In 

-8.520631 

0.100000 

-7.612530 

-7.607901 

0.200000 

-7.201041 

-7.266414 

0.300000 

-6.889720 

-6.979877 

0.400000 

-6.616004 

-6.715589 

0.500000 

-6.354047 

-6.454251 

0.600000 

-6.085607 

-6.177621 

0.700000 

-5.789374 

-5.857550 

0.800000 

-5.424880 

-5.419810 

0.900000 

-4.850462 

Table  5  95%  Confidence  Interval  Endpoints  (Model  2)  for  Test  Problem  1 


Left  Endpoint 
ln(XD) 

Probability 

Estimate 

Right  Endpoint 
ln{XD) 

-8.167631 

0.100000 

-7.660243 

-7.646688 

0.200000 

-7.233966 

-7.295645 

0.300000 

-6.917792 

-7.004956 

0.400000 

-6.642483 

-6.738876 

0.500000 

-6.380899 

-6.477293 

0.600000 

-6.114819 

-6.201983 

0.700000 

-5.824130 

-5.885810 

0.800000 

-5.473088 

-5.120757 

0.900000 

.4.449925 

Figure  4  95%  Confidence  Ellipse  for  Parameter  Estimates  a  and  />  in  Test  Problem  1 
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f 3  is  included  as  Figure  4.  Variances  have  also  been  computed  for  the  estimated  mean  and  standard 
deviation.  It  is  to  realize  that  these  parameters  are  calculated  for  the  transformed  explanatory 
variable,  XD  e  [-2,2].  The  estimated  mean  and  standard  deviation  are  computed  as  -0.820  and 
1.479,  respectively.  On  the  same  axis,  their  standard  deviations  are  respectively  estimated  to  be 
6.367  and  5.416. 
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3.2  Test  Problem  2:  Morbidity  of  Salmonella  T vphimuriuni 

This  test  problem  examines  the  morbidity  of  the  bacterium  Salmonella  Typhimurium  when 
administered  to  laboratory  mice.  A  reagent  containing  a  fraction  of  the  toxic  bacteria  was  prepared 
and  administered  to  seven  groups  of  mice. [21]  Each  group  consisted  of  five  mice.  In  this  case,  the 
number  of  survivors  is  reported  for  each  group  in  Table  6. 


Table  6  Salmonella  Typhimurium  Morbidity  Test  Data 


Dosage  (mg) 

Successes  out  of  5 

0.0625 

-3 

4 

0.125 

-2 

3 

0.25 

-1 

2 

0.5 

0 

0 

1.0 

1 

0 

2.0 

2 

0 

4.0 

3 

0 

The  calculations  entail  performing  a  fit  to  the  normal  curve  via  a  linear  argument  ( a  +  fix). 
Again,  Garwood’s  method  is  employed  as  in  Section  2.3.  From  a  casual  inspection  of  Table  6,  it 
is  evident  that  the  dosage  values,  after  natural  logarithmic  transformation,  conform  to  equation 

(3.1.1).  The  transformed  dosage  locus  X D  is  them  obtained  through  the  use  of  (3.1.2)  with  m  set 
to  three.  As  for  the  previous  test  problem,  the  starting  values  for  a  and  fi  are  kept  the  same  as  in 
equation  (3.1.3).  Garwood’s  algorithm  converges  to  the  final  values  below. 

a  =  1.51676;  fi  =  0.86344  (3.2.1) 

As  in  the  previous  case,  the  value  of  a  has  not  been  adjusted  by  the  addition  of  five  probits.  The 
addition  used  for  classical  analyses,  does  not  alter  the  results.  The  results  shown  here  agree  with 
the  archived  estimates  with  excellent  accuracy.  The  linear  fit  for  this  test  problem  is  presented  in 
Figure  5  and  is  compared  against  the  empirical  data.  The  attendant  probability  density  (PDF)  and 
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Figure  6  Probability  Distributions  for  Test  Problem  2 

cumulative  distribution  (CDF)  functions  are  presented  in  Figure  6.  The  endpoints  of  the  probability 
95%  confidence  interval  for  this  test  problem  are  provided  in  Table  7  for  Confidence  Model  1. 
Similar  information  is  provided  in  Table  8  for  Model  2,  the  standard  error  model.  Again  both  of 
these  intervals  are  calculated  for  the  95%  confidence  level.  The  notation  —  indicates  that  the 
predicted  endpoint  falls  out  of  the  domain.  Finally,  the  95%  confidence  ellipse  computed  for 
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Table  7  95%  Confidence  Interval  Endpoints  (Model  1)  for  Test  Problem  2 


Left  Endpoint 
ln(XD) 

Probability 

Estimate 

Right  Endpoint 

In 

-3.303797 

0.100000 

-2 .406570 

-2 .889236 

0.200000 

-2 .117651 

-2 . 605373 

0.300000 

-1 . 877926 

-2.370206 

0.400000 

-1 . 652407 

-2.155195 

0.500000 

-1 . 420311 

-1 . 943868 

0.600000 

-1 . 156728 

-1.720968 

0.700000 

-0.802384 

-1 .463119 

0.800000 

0 .442143 

-1 .107982 

0.900000 

— 

Table  8  95%  Confidence  Interval  Endpoints  (Model  2)  for  Test  Problem  2 


Left  Endpoint 
ln(X„) 

Probability 

Estimate 

Right  Endpoint 
ln(X0) 

-4.506577 

0.100000 

-2.588143 

-3 .110754 

0.200000 

-2.258334 

-2.749067 

0.300000 

-2 . 007947 

-2 .483716 

0.400000 

-1.784878 

-2.254286 

0.500000 

-1.567229 

-2.036637 

0.600000 

-1.337799 

-1 . 813568 

0.700000 

-1 . 072448 

-1.563181 

0.800000 

-0.710761 

-1.233372 

0.900000 

0.685063 

Figure  7  95%  Confidence  Ellipse  for  Parameter  Estimates  a  and  P  in  Test  Problem  2 
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estimates  a  and  f 3  is  shown  in  Figure  7.  The  mean  and  standard  deviation  for  the  fitted  distribution 
are  calculated  as  -1.756  and  1.158,  respectively,  while  the  corresponding  standard  deviations  for 
these  estimates  are  2.675  and  2.583.  Again,  these  values  are  computed  for  the  fully  transformed 

explanatory  variable  XD  e  [-3,3]. 
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3.3  Test  Problem  3:  Arsenic  Toxicity 


The  second  test  problem  examines  the  susceptibility  of  brine  shrimp  to  arsenic  in  liquid 
solution. [22]  As  the  concentration  of  arsenic  is  increased  through  a  geometric  progression,  we 
expect  fewer  shrimp  to  survive.  The  data  for  this  test  series  is  provided  in  Table  9.  The  true 

concentrations  for  arsenic  are  not  provided,  so  the  uniform  values  XD  are  used  as  listed. 


Table  9  Arsenic  Toxicity  Test  Data 


Solution 

Successes 

C 

-3 

8 

D 

-2 

8 

E 

-1 

6 

F 

0 

5 

G 

1 

5 

H 

2 

1 

I 

3 

0 

Figure  8  Normal  Curve  Linear  Fit  Form  for  Test  Problem  3 
In  this  case,  a  response  is  defined  as  a  brine  shrimp  surviving  immersion  in  the  solution. 

Our  starting  guess  for  the  estimation  routine  is  given  by  (3.2.1),  the  same  values  as  used  in  the 

preceding  test  problem.  The  converged  parameter  estimates  are 
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a  =  -0.434165;  J3  =  0.712834 


(3.3.1) 


These  values  constitute  an  excellent  match  for  the  archived  estimates,  and  the  final  linear  fit  is 
shown  in  Figure  8.  Since  no  actual  values  are  specified  for  the  solution  concentrations,  the  same 
dosage  levels  occurring  in  the  preceding  test  problem  are  used  here.  Probability  density  and 


Figure  9  Probability  Distributions  for  Test  Problem  3 


Table  10  95%  Confidence  Interval  Endpoints  (Model  1)  for  Test  Problem  3 


Left  Endpoint 

Probability 

Estimate 

Right  Endpoint 

-1 . 884088 

0.100000 

-1 . 021735 

-1.393996 

0.200000 

-0.660039 

-1 . 056911 

0.300000 

-0.370771 

-0.776747 

0.400000 

-0.106607 

-0.520028 

0.500000 

0.155914 

-0.267356 

0.600000 

0 . 438264 

-0.000706 

0.700000 

0 . 775710 

0.307552 

0.800000 

1.279043 

0.730878 

0.900000 

— 

cumulative  distribution  functions  are  easily  calculated  for  this  problem;  these  plots  are  shown  in 
Figure  9.  On  comparing  the  parameter  estimates  to  the  corresponding  values  in  Garwood,  one  may 
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Table  11  95%  Confidence  Interval  Endpoints  (Model  2)  for  Test  Problem  3 


Left  Endpoint 

Probability 

Estimate 

Right  Endpoint 

-2.240194 

0.100000 

-1 .166359 

-1.544010 

0.200000 

-0.767613 

-1 .161147 

0.300000 

-0 .467232 

-0.861901 

0.400000 

-0.201910 

-0.596125 

0.500000 

0.054181 

-0.340034 

0.600000 

0.319958 

-0 . 074711 

0.700000 

0.619203 

0.225669 

0.800000 

1.002066 

0 . 624415 

0.900000 

1 . 698250 

Figure  10  95%  Confidence  Ellipse  for  Parameter  Estimates  a  and  (3  in  Test  Problem  3 


note  that  the  probit  quantity  of  5  has  not  been  added  to  the  estimate  of  a  shown  in  (3.3.1). [11] 
With  this  adjustment  to  a  ,  the  agreement  is  excellent.  In  accordance  with  Confidence  Model  1, 
success  probability  confidence  intervals  are  provided  in  Table  10.  Corresponding  confidence 
intervals  computed  by  using  Model  2  are  shown  in  Table  11.  The  95%  confidence  ellipse  for  the 
estimates  of  a  and  /?  is  presented  in  Figure  10.  The  mean  and  standard  deviation  for  this  distribution 

are  0.609  and  1.402,  respectively.  These  parameters  are  computed  on  XD  the  transformed 
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explanatory  variable  axis  along  with  their  respective  standard  deviation  estimates  of  3.159  and 
3.221. 
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3.4  Test  Problem  4:  Dixon/Mood  Bruceton  First  Test  Series 


The  focus  of  this  test  problem  is  a  Bruceton  or  Up  and  Down  test.  As  described  earlier,  the 
test  consists  of  a  set  of  Bernoulli  trials;  each  trials  has  either  a  pass  or  fail  outcome.  This  problem 
is  taken  from  Dixon  and  Mood  and  is  comprised  of  60  trials  conducted  over  five  (dosage)  levels 
in  the  explanatory  variable.  The  test  data  can  be  formulated  in  terms  of  five  Binomial  trials.  The 
data  in  this  form  is  provided  in  Table  12. 


Table  12  Bruceton  Test  Data  for  Test  Problem  4 


Dosage  Level 

Number  of  Trials 

Successes 

2.0 

1 

1 

1.7 

10 

10 

1.4 

27 

18 

1.1 

20 

2 

0.8 

2 

0 

For  this  test  problem,  the  algorithm  described  in  Section  2.6  is  employed  to  estimate  the  mean  and 
standard  deviation  for  the  normal  probability  distribution.  Initial  values  for  the  mean  /u  and 
standard  deviation  a  are  set  as  follows: 

ju0  =1.4  cr0  =0.18  (3.4.1) 


Figure  1 1  Probability  Distributions  for  Test  Problem  4 
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The  iterative  scheme  is  convergent  rendering  the  final  estimates 


H  =  1.3216;  a- =  0.1667  (3.4.2) 

With  these  values  in  mind,  the  probability  distributions  for  this  test  problem  are  shown  in  Figure 
11.  95%  confidence  intervals  for  the  success  probability  according  to  Confidence  Model  1  are 
documented  in  Table  13.  Similar  information  for  Confidence  Model  2  is  contained  in  Table  14. 
The  confidence  ellipse  for  parameters  //  and  u  is  shown  in  Figure  12.  The  standard  deviations 
Table  13  95%  Confidence  Interval  Endpoints  (Model  1)  for  Test  Problem  4 


Left  Endpoint 

Probability 

Estimate 

Right  Endpoint 

1 . 0467469 

0.100000 

1 .  1897100 

1 .  1305151 

0.200000 

1.2520503 

1 . 1881752 

0.300000 

1.3016885 

1 . 2361231 

0.400000 

1 . 3468707 

1 .2800741 

0.500000 

1.3916125 

1 . 3233422 

0.600000 

1 . 4394931 

1.3690090 

0.700000 

1.4961876 

1 . 4217988 

0.800000 

1 . 5784431 

1 .4942660 

0.900000 

— 

Table  14  95%  Confidence  Interval  Endpoints  (Model  2)  for  Test  Problem  4 


Left  Endpoint 

Probability 

Estimate 

Right  Endpoint 

0.9921593 

0.100000 

1 .1663933 

1 . 1067824 

0.200000 

1 . 2347979 

1 .  1715661 

0.300000 

1 .2862759 

1 . 2225015 

0.400000 

1.3316982 

1 .2678664 

0.500000 

1 . 3754852 

1.3116534 

0.600000 

1 . 4208501 

1 . 3570757 

0.700000 

1 . 4717855 

1.4085537 

0.800000 

1.5365692 

1.4769583 

0.900000 

1 . 6511923 
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Figure  12  95%  Confidence  Ellipse  for  Parameter  Estimates  p  and  a  in  Test  Problem  4 

for  n  and  a  are  computed  from  the  asymptotic  variance-covariance  as  0.1664  and  0.1505, 
respectively. 
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3.5  Test  Problem  5:  Dixon/Mood  Bruceton  Second  Test  Series 

This  section  of  the  report  discusses  the  second  Mood/Dixon  Bruceton  test  problem.  As  for 
the  preceding  case,  this  test  entails  60  trials  taken  over  four  levels  of  the  explanatory  variable.  The 
test  data,  compiled  as  a  series  of  four  binomial  trials,  is  provided  in  Table  15. 


Table  15  Bruceton  Test  Data  for  Test  Problem  5 


Dosage  Level 

Number  of  Trials 

Successes 

1.9 

3 

3 

1.5 

29 

27 

1.9 

27 

1 

0.1 

1 

0 

Figure  13  Probability  Distributions  for  Test  Problem  5 


The  initial  estimates  for  /u  and  a  are  the  same  as  in  (3.4.1).  The  iterative  scheme  converges  to  the 
final  values  of 

fi  =  1.22775;  <r  =  0.18347  (3.5.1) 

Accordingly  the  probability  density  and  cumulative  distribution  functions  for  the  normal 
distribution  fit  to  the  data  is  shown  in  Figure  13.  Confidence  intervals  for  the  success  probability 
are  listed  in  Table  16  for  Confidence  Model  1  and  in  Table  17  for  Confidence  Model  2.  The 
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Table  16  95%  Confidence  Interval  Endpoints  (Model  1)  for  Test  Problem  5 


Left  Endpoint 

Probability 

Estimate 

Right  Endpoint 

0.9252993 

0.100000 

1 . 0825744 

1 . 0174535 

0.200000 

1 .  1511556 

1 . 0808859 

0.300000 

1.2057630 

1 . 1336339 

0.400000 

1 . 2554684 

1 .1819848 

0.500000 

1.3046893 

1 . 2295845 

0.600000 

1.3573632 

1 . 2798229 

0.700000 

1 .  4197333 

1.3378975 

0.800000 

1 . 5102234 

1 .4176194 

0.900000 

— 

Table  17  95%  Confidence  Interval  Endpoints  (Model  2)  for  Test  Problem  5 


Left  Endpoint 

Probability 

Estimate 

Right  Endpoint 

0.8652470 

0.100000 

1 . 0569234 

0.9913450 

0.200000 

1 . 1321760 

1 . 0626141 

0.300000 

1 . 1888074 

1 .1186486 

0.400000 

1.2387769 

1 .1685550 

0.500000 

1 .2869475 

1 .2167255 

0.600000 

1.3368538 

1.2666950 

0.700000 

1 .3928884 

1 . 3233265 

0.800000 

1 .4641575 

1.3985790 

0.900000 

1.5902555 

1.15  1.2  jj  1.25  1.3 

Figure  14  95%  Confidence  Ellipse  Computed  for  Parameter  Estimates  /./  and  n  in  Test  Problem  5 
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confidence  ellipse  for  the  estimated  values  of  //  and  a  is  shown  in  Figure  14.  Standard  deviation 


approximations  have  been  calculated  for  the  estimated  p  and  a.  Respectively,  these  values  are 
0.27238  and  0.16754. 
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3.6  Test  Problem  6:  Anti-Pneumococcus  Serum 


In  the  structure  of  this  report,  this  test  problem  is  new.  This  report’s  predecessor  [16]  did 
not  discuss  the  logistic  probability  distribution.  Yet,  this  statistical  distribution  is  widely  applied, 
so  it  has  be  included  in  this  revised  report.  This  test  problem  is  taken  from  Collett  [4]  and  represents 
a  means  of  testing  antibiotic  dosage  efficacy.  In  this  clinical  test  series,  five  antibiotic  dosage  levels 
are  selected.  For  each  level,  40  mice  are  infected  with  pneumococcal  bacterial  and  then  injected 
with  the  dosage  for  the  level.  Then  the  number  of  mice  that  do  not  survive  after  a  seven  days 
waiting  period  at  that  level  are  recorded.  The  binomial  data  for  this  test  series  is  presented  in  Table 
18.  A  success  in  this  case  is  interpreted  as  a  non-survivor.  The  logistic  distribution  is  used  to  fit 


Table  18  Logistic  Test  Data  for  Test  Problem  6 


Dosage  Level 

Number  of  Trials 

Successes 

0.0028 

40 

35 

0.0056 

40 

21 

0.0112 

40 

9 

0.0225 

40 

6 

0.0450 

40 

1 

Figure  15  Linear  Logistics  Fit  Form  for  Test  Problem  6 
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Figure  16  Probability  Distributions  for  Test  Problem  6 


Table  19  95%  Confidence  Interval  Endpoints  (Model  1)  for  Test  Problem  6 


Left  Endpoint 
ln(Xfl) 

Probability 

Estimate 

Right  Endpoint 
ln(Xfl) 

-3.4005 

0.100000 

-4 . 3251 

-3.9572 

0.200000 

—4 . 6660 

-4 .2990 

0.300000 

-4 . 9281 

-4.5654 

0.400000 

-5.1651 

-4 .8004 

0.500000 

-5.4031 

-5.0276 

0.600000 

-5.6671 

-5.2675 

0.700000 

-6.0009 

-5.5514 

0.800000 

-6.5509 

-5.9659 

0.900000 

— 

this  data  with  the  natural  logarithm  of  the  dosage  as  an  explanatory  variable.  Moreover,  the  linear 
fitting  form  (2.9.1)  is  employed.  At  the  start  of  the  iterative  scheme,  the  following  values  are 
assumed  for  a  and  ft. 


aQ  =  -9.19;  /?0=-1.83 


(3.6.1) 
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Table  20  95%  Confidence  Interval  Endpoints  (Model  2)  for  Test  Problem  6 


Left  Endpoint 
ln(XD) 

Probability 

Estimate 

Right  Endpoint 
ln(XD) 

-2 . 9558 

0.100000 

-4 .1853 

-3.7993 

0.200000 

-4.5689 

-4 .1985 

0.300000 

-4 .8433 

-4 .4875 

0.400000 

-5.0815 

-4 .7328 

0.500000 

-5.3123 

-4 . 9637 

0.600000 

-5.5576 

-5.2019 

0.700000 

-5.8466 

-5.4763 

0.800000 

-6.2458 

-5.8598 

0.900000 

-7.0893 

Figure  17  95%  Confidence  Ellipse  for  Parameter  Estimates  a  and  [>  in  Test  Problem  6 

The  iterative  scheme  converges  to  the  final  parameter  estimates  of 

a  =  -9.1893;  £  =  -1.8296  (3.6.2) 

The  attendant  curve  fit  is  shown  in  Figure  15  while  the  probability  density  and  cumulative 
distribution  functions  are  shown  in  Figure  16.  Based  upon  Confidence  Model  1,  the  confidence 
interval  endpoints  for  success  probability  are  listed  in  Table  19.  Similar  information  for 
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Confidence  Model  2  is  provided  in  Table  20.  The  confidence  ellipse  for  this  test  problem  is 
presented  in  Figure  17.  The  mean  and  standard  deviation  for  this  distribution  are  easily  computed 
from  the  linear  form  (2.3.2)  as 

//  = -5.02256;  a  =  -0.54656  (3.6.3) 

based  upon  the  logarithm  of  the  explanatory  variable.  From  the  asymptotic  variance-covariance 
matrix,  the  standard  deviations  for  the  parameter  estimates  (3.6.3)  are  computed  as  9.69016  and 
13.40990,  respectively. 
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3.7  Test  Problem  7:  Aircraft  Fastener  Failure 


This  study  involves  the  failure  of  aircraft  fasteners  under  compressive  loading.  Series  of 
tests  were  performed  on  collections  of  fasteners  for  ten  different  compressive  loads  measured  in 
psi.  For  each  fastener  collection  (or  lot),  a  certain  number  of  fasteners  fail  resulting  in  a  binomial 
series  with  load  as  the  explanatory  variable.  Table  21  contains  the  data  for  the  series  of  fastener 
tests.  For  this  test  case,  the  load  is  directly  used  as  the  explanatory  variable  without  transformation. 
As  in  the  previous  test  problem,  the  logistic  distribution  is  applied.  The  resulting  linear  fit  is 


Table  21  Logistic  Test  Data  for  Test  Problem  7 


Load  (PSI) 

Lot  Size 

Number  Failing 

2500 

50 

10 

2700 

70 

17 

2900 

100 

30 

3100 

60 

21 

3300 

40 

18 

3500 

85 

43 

3700 

90 

54 

3900 

50 

33 

4100 

80 

60 

4300 

65 

51 

Figure  18  Linear  Logistics  Fit  Form  for  Test  Problem  7 
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Figure  19  Probability  Distributions  for  Test  Problem  7 


shown  in  Figure  18.  The  probability  density  and  cumulative  distribution  functions  for  the  logistic 
distribution  are  presented  in  Figure  19.  In  accordance  with  Confidence  Model  1,  the  success  ratio 
confidence  interval  endpoints  are  provided  in  Table  22.  Corresponding  endpoints  for  Model  2  are 
provided  in  Table  23.  The  confidence  ellipse  for  distribution  parameters  a  and  fi  is  shown  in 
Table  22  95%  Confidence  Interval  Endpoints  (Model  1)  for  Test  Problem  7 


Left  Endpoint 

Probability 

Estimate 

Right  Endpoint 

1745 .1869 

0.100000 

2341.7176 

2343.7843 

0.200000 

2794 .1246 

2721.8920 

0.300000 

3116.8377 

3022 .1416 

0.400000 

3393.5079 

3290 . 9217 

0.500000 

3657.2709 

3553.8612 

0.600000 

3931.5218 

3834 .3142 

0.700000 

4244.8644 

4168 . 6696 

0.800000 

4655.5562 

4657 . 8710 

0.900000 

5387 .4483 

Figure  20.  As  is  evident,  the  fitting  procedure  functions  quite  well  rendering  the  parameter 
estimates  of 
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■3 


(3.7.1) 


a  =  -5.33971  and  /?  =  1.54843  x  10“ 


Table  23  95%  Confidence  Interval  Endpoints  (Model  2)  for  Test  Problem  7 


Left  Endpoint 

Probability 

Estimate 

Right  Endpoint 

1650.0926 

0.100000 

2285.9406 

2298 .4357 

0.200000 

2758 .7448 

2690 .1111 

0.300000 

3087 . 8241 

2995.9883 

0.400000 

3366.5277 

3267.1356 

0.500000 

3629.7845 

3530.3924 

0.600000 

3900.9318 

3809.0960 

0.700000 

4206.8090 

4138 .1752 

0.800000 

4598 .4844 

4610 . 9794 

0.900000 

5246.8275 

Figure  2  95%  Confidence  Ellipse  for  Parameter  Estimates  a  and  /?  in  Test  Problem  7 


The  corresponding  mean  and  standard  deviation  for  this  distribution  are  computed  as 

fi  =  3448.460  and  cr  =  645.813  (3.7.2) 

A  variance  can  be  calculated  for  these  two  parameters,  i.e.,  1.86844  xlO2  and  1.52363x10  2 , 
respectively. 
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4  CONCLUSIONS 


In  this  report,  we  have  presented  a  discussion  of  sensitivity  (or  Go/No-Go)  testing  from 
basic  principles.  The  prototypical  set  of  sensitivity  experiments  is  based  upon  the  Bruceton  test 
procedure,  or  the  “Up  and  Down  Method”,  used  to  estimate  the  statistical  parameters  associated 
with  a  series  of  binary  trials.  The  outcome  of  each  binary  trial  is  a  random  variable,  and  these 
random  variables  are  independent  and  identically  distributed.  For  the  “Probit”  method,  the  mean 
and  standard  deviation  are  estimated  by  fitting  data  to  a  normal  distribution.  For  pharmaceutical 
applications,  the  so-called  “dosage-mortality  curve”  is  produced  in  the  form  of  a  linear  fit.  The 
behavior  of  the  normal  curve  is  then  used  to  determine  the  probability  of  success  associated  with 
an  explanatory  variable  of  choice.  For  the  Logit  method,  the  logistic  distribution  is  used  since  it  is 
a  bit  more  conservative  near  the  tails,  yet  it  looks  a  lot  like  the  normal  distribution.  The  set-up  and 
execution  of  the  Bruceton  test  procedure  have  been  discussed  for  an  example  scenario.  The 
importance  of  locating  the  mean  is  stressed  as  well  is  ensuring  that  extraneous  data  far  removed 
from  the  mean  is  excluded  from  the  analysis.  Moreover,  issues  surrounding  the  efficiency  of  the 
data  have  been  discussed.  That  is  to  say,  we  must  obtain  a  sufficiently  large  sample  in  order  for 
our  analysis  techniques  and  assumptions  to  apply.  Algorithms  for  data  analysis  have  been 
discussed  in  detail,  particularly  on  the  topics  of  Probit  analysis,  Logit  analysis  and  confidence 
estimation.  We  have  shown  the  set  of  equations  required  to  fit  go-no  go  data  to  the  normal 
distribution  and  to  the  logistic  distribution.  These  methods  have  proven  capable  and  they  readily 
address  binomial  trials  on  nonuniformly  separated  levels.  The  fundamental  solution  scheme 
discussed  here  is  based  upon  Garwood’s  method,  yet  we  have  presented  an  alternative  to  this 
method  that  can  easily  incorporate  both  binary  and  binomial  trials.  The  algorithm  structure  works 
for  both  the  Probit  and  Logit  distributions.  The  alternate  algorithm  is  stable  and  remains  well- 


Distribution  A 
83 


defined  in  every  case.  From  the  standpoint  of  sampling,  we  have  discussed  interval  estimation  for 
sensitivity  tests.  Algorithms  have  been  suggested  for  the  qualified  estimation  of  probability 
confidence  intervals  as  well  as  those  for  the  mean  and  variance.  We  have  also  derived  equations 
that  allow  the  variance  and  covariance  to  be  estimated  for  the  normal  mean  and  variance  produced 
by  our  fitting  procedure.  Seven  classic  test  problems  have  been  chosen  from  engineering  and 
biological  science  to  serve  as  practical  examples.  These  problems  have  been  solved  by  directly 
coding  the  algorithms  presented  in  this  report.  In  each  case,  our  results  have  achieved  excellent 
agreement  with  archival  solutions.  Confidence  regions  have  been  computed  for  these  cases.  We 
have  also  discussed  shortfalls  that  exist  within  our  process  for  estimating  intervals  and  for 
estimating  elements  in  the  asymptotic  variance-covariance  matrix. 

Sensitivity  testing  is  very  important,  not  only  to  the  government  sector,  but  also  to  the 
commercial  sector,  namely  the  pharmaceutical  industry.  The  efficacies  and  effective  dosages  for 
medications  are  largely  determined  through  this  type  of  testing  and  analysis.  In  most  cases,  the 
statistical  distributions  produced  by  the  analysis  are  used  to  make  important  management 
decisions.  The  investment  in  the  development  of  a  single  medication  can  cost  hundreds  of  millions 
of  dollars.  The  cost  of  qualifying  explosive  systems  is  also  high.  For  these  reasons,  it  is  important 
that  we  attain  a  good  understanding  of  the  probability  and  statistics.  It  is  very  easy  to  become  mired 
in  either  of  these  fields.  These  divisions  of  mathematics  undergird  sensitivity  testing  and  analysis. 
System  reliability  is  of  paramount  importance  for  protecting  both  the  investment  of  funding  and 
human  life.  Failing  to  accurately  estimate  the  reliability  of  a  modem  drug  or  engineering  system 
can  have  grave  consequences. 
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